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I.  INTRODUCTION  AND  SUMMARY 


1.1  BACKGROUND 

To  successfully  monitor  a  test  ban  treaty  prohibiting  or 
limiting  underground  nuclear  exolosions,  it  is  necessary  to  under¬ 
stand  the  seismic  signatures  of  these  events.  An  important  part  of 
the  research  effort  to  imorove  this  understanding  has  been  the 
development  and  apolication  of  deterministic  methods  to  compute  the 
seismic  wave  signatures  of  nuclear  explosions.  In  nearly  all  of 
this  work  the  source  is  assumed  to  be  spherically  symmetric  and  the 
earth  is  assumed  to  be  plane-layered.  Theoretical  seismograms  can 
then  be  computed  with  widely  available  methods  and  .omoared  to 
observations.  This  procedure  has  been  quite  successful  and  most  of 
the  important  controlling  parameters  have  been  identified  and  their 
effects  have  been  quantified  to  some  degree.  However,  many 
important  issues  remain  unresolved. 

Most  of  the  important  questions  that  remain  regarding  the 
generation  of  seismic  waves  by  underground  explosions  are  associated 
with  multi-dimensional  effects.  For  example,  there  is  not  yet  a 
clear  understanding  of  the  effect  on  the  seismic  radiation  of  nearby 
interfaces,  the  free  surface  (allowing  spallation),  the  overburden 
pressure,  nonhydrostatic  prestress,  and  zones  of  weakness  in  the 
near  source  environment.  The  one-dimensional  calculations  now 
incorporate  detailed  constitutive  models  that  include  realistic 
models  for  pore  collaose,  effective  stress,  yielding,  and  cracking 
due  to  shear  and  tension  failure.  These  constitutive  models,  have 
been  generalized  to  two-dimensions  for  the  axisymmetric  finite 
difference  source  calculations  analyzed  ’o  this  nsoort.  While  they 
do  not  include  all  the  multi-dimensional  effects  we  have  listed, 
they  do  include  some  of  the  most  important,  and  so  reDresent  a 
significant  step  forward  in  the  development  of  realistic  theoretical 
simulations  of  the  seismic  waves  from  underground  explosions. 
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1.2  AXISYMMETRIC  CALCULATIONS  OF  NUCLEAR  EXPLOSIONS 


In  th^s  report  we  present  a  detailed  analysis  of  the  seismic 
waves  from  eleven  two-dimensional  finite  difference  calculations  of 
underground  nuclear  explosions  in  granite.  Seven  of  these  calcu¬ 
lations  were  done  by  J.  Trulio  and  N.  Perl  of  Applied  Theory,  Inc. 
(ATI)  and  four  by  N.  Rimer  and  J.  T.  Cherry  of  Systems,  Science  and 
Software  (S-Cubed).  They  share  the  same  axi  symmetric  geometry, 
though  the  details  are  different. 

The  two-dimensional,  axisymmetric  exDlosion  calculations 
include  the  presence  of  a  free  surface,  the  dependence  of  overburden 
pressure  on  depth  and,  for  the  S-Cubed  calculations,  some  dependence 
of  material  properties  on  depth.  These  are  probably  the  most 
important  higher  order  corrections  to  the  one-dimensional  models. 
The  depth-dependence  of  overburden  pressure  is  a  property  of  all 
test  sites.  Further,  our  understanding  of  aeolcgic  structures  is 
normally  based  on  plane-layered  models.  Therefore,  axisymmetry  is  a 
natural  geometry,  and  specification  of  the  geometry  and  material 
properties  in  two-dimensions  is  more  straightforward  than 
characterizing  the  entire  near-source  environment  by  a  homogeneous 
material . 

The  most  important  two-dimensional  effects  are  associated  with 
the  nonlinear  interaction  of  the  stress  waves  with  the  free  sur¬ 
face.  Surface  spallation  is  an  obvious,  even  dominant,  phenomenon 
observed  in  the  near-field,  yet  it  has  never  been  included  in 

seismic  wave  oropagation  studies  in  a  very  satisfactory  way. 

Perhaps  more  imoortant,  there  is  amole  evidence  that  the  free 

surface  ohase  oP  is  more  complex  than  predicted  by  elastic  theory 

with  soherically  symmetric  sources,  but,  again,  this  remains  mostly 
in  the  realm  of  speculation.  With  the  axisymmetric  calculations,  we 
are  able  to  study  these  important  effects. 

The  objectives  of  the  parameter  variations  in  the  ATI  and 
S-Cubed  suites  of  calculations  are  somewhat  different.  The  ATI 
calculations  were  all  done  at  the  same  yield,  150  kt,  in  a 
hypothetical  granite  halfspace.  Only  the  depth  was  varied,  from  159 
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to  1000  meters.  The  philosophy  of  the  S-Cubed  calculations  was  to 
begin  by  modeling  a  specific  event,  PILEDRIVER,  at  the  Nevada  Test 
Site.  The  constitutive  model  and  source  geology  (three  layers)  were 
chosen  for  this  event.  The  computed  and  observed  ground  motions 
were  compared  at  some  twenty-five  near-field  gauge  locations.  The 
comparison  was  quite  good,  except  that  the  calculation  overpredicted 
the  amount  of  cracking  (spallation)  within  300  meters  of  ground 
zero.  This  turns  out  to  be  important  when  comparing  synthetic  and 
observed  far-field  seismograms. 

The  other  three  S-Cubed  source  calculations  were  the  same, 
except  the  depth  and  yield  were  varied  from  the  values  appropriate 
for  PILEDRIVER.  The  four  calculations  were  as  follows: 

Depth  (Meters)  Yield  (Kt) 


463 

1000 

1000 

400 


60  (PILEDRIVER) 
150 
20 
20 


1.3  OUTLINE  OF  THE  ANALYSIS 

The  analyses  of  the  seven  ATI  and  four  S-Cubed  calculations 
are  described  separately,  but  follow  parallel  lines.  First,  we 
describe  the  near-field  ground  motions  predicted  by  the  calculations 
in  some  detail.  Then  we  present  theoretical  seismograms  for  these 
calculations  and  analyze  their  implications.  This  outline  is 
aoparent  in  the  section  headings  listed  below: 


Section  II:  ATI  Granite  Calculations 

Section  III:  Ms  and  m^,  Estimates  for  ATI 
Granite  Calculations 


Section  IV:  S-Cubed  Granite  Calculations 

Section  V:  N!s  and  m^  Estimates  for  S-Cubed 

Granite  Calculations 
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Section  VI:  Comparison  of  Predicted  and  Observed 

Seismograms  for  PILEDRIVER 

Section  VII:  Comparative  Analysis  of  ATI  and 
S-Cubed  Granite  Calculations. 

A  key  step  in  this  study  is  the  linkage  between  the  near¬ 
field  ground  motions  computed  by  the  finite  difference  programs  and 
the  analytical  techniques  used  to  propagate  seismic  waves  in 
realistic  earth  models.  A  rigorous  method  for  accomplishing  this 
linkage  is  described  in  Appendix  A,  "Synthetic  Seismograms  from 
Complex  Source  Calculations."  The  details  and  examples  presented  in 
this  appendix  are  mostly  for  the  calculation  of  the  normal  mode 
(Rayleigh  waves)  displacements  in  a  plane-layered  earth  model,  but 
ray  theory  methods  for  propagating  body  waves  can  be  used  within  the 
same  theoretical  framework. 

To  compute  synthetic  seismograms  at  large  distances  for  the 
finite  difference  source  calculations,  it  is  necessary  to  monitor 
the  tractions  and  displacements  on  some  (hypothetical)  surface  which 
entirely  encloses  the  region  of  nonlinear  material  response.  Our 
discussion  of  the  source  calculations  in  Sections  II  and  IV  is 
concerned  with  the  characteristic?  of  the  monitored  ground  motions 
on  this  elastic  surface. 

An  important  constraint  on  the  numerical  results  is  the  re¬ 
quirement,  based  on  conservation  of  momentum,  that  the  total  down¬ 
ward  force  and  impulse  (on  the  surface  enclosing  the  source  region) 
vanish  at  late  time.  It  turns  out  (Appendix  A)  that  the  inevitable 
numerical  errors  that  cause  deviation  from  tnis  requirement  can 
dominate  the  solution  at  long  periods.  Therefore,  a  major  theme  in 
Sections  II  and  IV  is  the  application  of  a  "correction"  to  adjust 
the  computed  vertical  force  and  impulse  to  zero  at  the  last  time 
step. 

The  synthetic  seismogram  results  are  Dresented  in  Sections  III 
and  V.  The  vertical  force  "correction"  is  also  an  issue  in  these 
sections,  because  it  would  be  unsatisfactory  for  an  ad_  hoc  cor¬ 
rection,  which  is  what  we  apply,  to  dominate  the  answer.  Our 
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conclusion  is  that  the  vertical  force  terms,  after  correction,  do 
not  play  an  important  role.  All  of  the  solution  checks  we  are  able 
to  make  indicate  that  the  body  and  surface  wave  synthetic 
seismograms  accurately  represent  the  seismic  waves  generated  by 
these  theoretical  sources. 

The  synthetic  seismogram  analyses  in  Sections  III  and  V  are 
focussed  on  comparing  the  seismic  waves  from  the  two-dimensional 
sources  with  those  from  analogous  one-dimensional,  reduced  displace¬ 
ment  potential  (ROP)  sources.  For  the  S-Cubed  calculations,  we  are 
able  to  compare  to  the  RDP  source  computed  with  the  same  constitu¬ 
tive  model  in  spherical  symmetry.  No  such  one-dimensional  calcu¬ 
lation  was  available  for  the  ATI  granite,  so  we  compare  to  the  ROP 
predicted  by  the  semi -empirical  model  of  Mueller  and  Murphy  (1971). 

Most  of  the  calculations  are  for  hypothetical  events,  so 
direct  comparison  with  observed  seismograms  is  only  possible  for  the 

PILEDRIVER  calculation  done  by  S-Cubed.  This  comparison  is  made  in 

Section  VI.  Finally,  in  Section  VII,  the  results  of  all  eleven 

calculations  are  plotted  together  for  direct  comoarison. 

1.4  CONCLUSIONS 

This  study  is  primarily  an  investigation  of  the  influence  of 
burial  depth  on  the  seismic  signals  from  underground  explosions. 
For  the  most  part,  the  results  are  interesting,  but  not  terribly 
exciting,  because  they  are  pretty  much  in  accord  with  our  expec¬ 
tations.  The  key  exception  is  the  surface  wave  results  for  the 
S-Cubed  calculations. 

"ne  most  important  result  of  this  study  is  that  the  S-Cubed 
calculations  show  M  to  be  a  strong  function  of  depth.  The 

shallow  S-Cubed  calculations  have  surface  wave  amplitudes  that  are  a 
factor  of  two  or  three  larger  than  those  from  a  comparable  one-di¬ 
mensional  source  calculation.  This  effect  is  probably  exaggerated 
because  the  free  surface  interaction  effects  are  too  large  in  the 
S-Cubed  calculations,  but  would  remain  important  even  if  the  near 
surface  material  were  strengthened.  On  the  other  hand,  the  ATI 
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calculations  show  no  strong  dependence  of  M  on  depth,  even  thouah 
the  shallow  ATI  sources  cratered.  U'e  do  not  know  whv  the  two  sets 
of  calculations  give  such  different  results. 

The  important  conclusions  are  listed  at  the  end  of  several 
sections  of  the  report.  Summaries  for  the  ATI  and  S-Cubed  calcu¬ 
lations  are  given  separately  in  Sections  3.7  and  5.10.  Cur  conclu¬ 
sions  about  the  comparison  of  synthetic  and  observed  seismograms  are 
listed  in  Section  6. A,  and  all  of  Section  VII  should  be  read  as  a 
summary. 

The  main  results  are  listed  below: 

ATI  Calculations 

•  Neither  mb  nor  are  strongly  dependent  on 
depth.  The  most  effect  was  on  M  at  shallow 
depths. 

•  Compared  to  the  Mueller /Murphy  RDP  source,  the 
depth  dependence  in  both  amolitude  and  corner 
frequency  is  less  for  the  ATI  sources. 

•  The  pP  phase  appears  to  be  smaller  than 
expected  from  elastic  theory. 

S-Cubed  Calculations 

•  Two-dimensional  effects  are  not  very  important 

for  the  two  deep  explosions  (20  kt  and  150  kt 
at  1000  meters).  Both  and  are  little 

different  from  the  values  estimated  from  an  RDP 
source  computed  with  the  same  constitutive 
model  at  the  same  depth.*  This  is  true  even 
though  considerable  cracking  and  spallation 
occur  in  the  two-dimensional  calculations. 

*  These  results  are  essentially  tests  of  the  entire  computational 
procedure.  It  is  gratifying  that  these  two  very  different  and 
complex  procedures  arrive  at  the  same  results. 
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•  For  body  waves  the  first  arriving  P  wave  is 
essentially  the  same  for  one-  and  two-di¬ 
mensional  sources.* 

•  The  two-dimensional  effects  enhance  the  surface 
wave  amplitudes  for  the  shallow  events  (60  kt 
at  463  meters  and  20  kt  at  400  meters)  by  a 
factor  of  two  or  three.  We  must  qualify  this 
by  pointing  out  that  these  shallow  calculations 
have  very  strong  surface  interaction  effects. 

Comparison  with  PHEORIVER  data  indicates  that 
the  free  surface  interaction  is  overpredicted, 
at  least  for  that  event.  Less  free  surface 
interaction  would  presumably  give  less 
enhancement  of  the  surface  wave  amplitudes. 

•  For  the  shallow  events,  the  m£  is  dif¬ 
ferent  than  predicted  with  an  RDP  source, 
though  by  less  than  0.2  units. 

•  Analysis  of  the  spectra  show  that  in  no  case  is 
pP  a  spectral  shadow  of  ?,  as  it  is  for  an  RDP 
source  and  elastic  propagation. 

•  Phases  that  seem  to  be  associated  with  spall 
closure  can  be  seen  on  the  shallow  source  body 
wave  records.  However,  they  are  not  easily 
associated  with  identifiable  crack  closure 
patterns  in  the  source  calculation. 

The  comoarison  of  observed  and  calculated  body  and  surface 
waves  for  PILcDRI’/ER  in  Section  '/I  leads  to  the  conclusion  that  the 
two  are  in  rather  good  agreement.  This  conclusion  is  subject  to  the 
qualifications  one  often  faces  in  this  kind  of  comparison.  For  the 
surface  waves,  it  is  the  need  to  account  for  the  non-axi symmetric 

*  See  footnote  on  orevious  page. 
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comoonent,  usually  attributed  to  tectonic  stress  release,  Medina  a 
recent  estimate  for  this  component  by  Rivers  and  von  Seggern  (1979) 
to  our  solution,  we  get  good  agreement  with  the  data. 

For  body  waves  the  comparison  is  complicated  by  the  apparent 
presence  of  strong  azimuthal  effects  in  the  radiated  short  period 

energy  (Hadley  and  Hart,  1979).  However,  our  conclusion  is  that  the 
computed  PILEQRIVER  source  (in  one-  or  two-dimensions)  has  about  the 
right  direct  P  amplitude.  However,  the  two-dimensional  source 
calculation  appears  to  include  too  much  non-linear  interaction  with 
the  free  surface.  The  later  portion  of  the  P  waveform  does  not 

match  the  data,  apparently  because  pP  is  too  greatly  suooressed  and 
because  the  seismic  energy  from  spall  closure  is  too  large  or  is 
timed  incorrectly.  This  "overprediction"  of  surface  interaction 
effects  is  expected  since  comparison  of  theoretical  and  observed 

near-field  motions  and  plots  of  the  cracking  near  the  source 
indicate  that  there  was  too  much  spallation  in  the  calculation. 

All  the  theoretical  and  m^  values  are  plotted  together 

in  Section  VII.  They  are  shown  versus  source  depth  and  versus 
yield,  including  the  observed  values  for  HARDHAT,  SHOAL  and  °ILE- 
DRI'/ER.  The  most  dramatic  difference  between  the  ATI  and  S-Cubed 
calculations  is  the  strong  dependence  of  M  on  depth  predicted  by 
S-Cubed,  but  not  by  ATI  sources.  This  must  be  a  reflection  of  the 
different  constitutive  models  used. 

The  constitutive  models  used  by  S-Cubed  (in  spherically 
symmetric  source  calculations)  lead  to  RDP  source  functions  that  are 
strongly  oeaked,  with  the  value  near  1  uz  a  factor  of  five  or  more 
larger  than  the  value  at  long  periods.  rne  oeaking  is  due  to  the 
incorporation  of  an  effective  stress  law  ana  the  choice  of 
unconfined  compressive  strength  (0.75  kbar),  based  on  laboratory 
data  for  fractured  granite  and  results  of  comparison  with  near-field 
ground  motion  observations  (Rimer,  personal  communication). 

While  we  do  not  have  RD°  source  functions  for  the  ATI  granite, 
comparison  of  M  and  m^  for  the  ATI  two-dimensional  calculations 
indicates  that  the  RCP  peaking  is  probably  less  than  a  factor  of 
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two.  Due  primarily  to  this  difference,  the  m^  for  the  S-Cubed 
calculations  is  about  0.5  units  higher  than  that  for  ATI  calcula¬ 
tions  of  the  same  yield.  The  ATI  M  values  fall  between  those  for 

s 

the  shallow  and  deep  S-Cubed  calculations. 
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II.  ATI  GRANITE  CALCULATIONS 


2.1  INTRODUCTION 

We  will  subsequently  be  analyzing  the  output  of  seven  source 
calculations  done  by  J.  Trulio  and  N.  Perl  of  Applied  Theory,  Inc. 
The  purpose  of  this  analysis  is  to  determine  the  body  and  surface 
wave  magnitudes  (m^  and  M  )  associated  with  the  simulated 

seismic  events,  and  this  is  done  with  the  method  described  in 
ADoendix  A. 

Our  intention  in  this  section  is  to  describe  the 
characteristics  of  the  source  calculations  and  the  comouted 
near-source  ground  motions.  An  important  constraint  on  the 
numerical  results  is  the  requirement,  based  on  conservation  of 
linear  momentum,  that  the  total  downward  force  and  impulse  vanish  at 
late-time.  Due  to  numerical  errors  that  are  inevitable  in  source 
calculations  of  this  kind,  this  requirement  is  not  satisfied 

exactly.  The  data  are  therefore  "corrected"  by  imposing  an 
additional  time  step  that  causes  the  force  and  impulse  to  vanish. 
This  correction  is  described  in  Section  2. a.  It  is  the  "corrected" 
data  that  are  used  for  the  synthetic  seismogram  calculations 
described  in  Section  III, 


2.2  SOURCE  DESCRIPTION 

The  seven  ATI  calculations  are  for  150  KT  explosions  in  a 
granite  half space  characterized  by 


p  wave  velocity:  a  =  a. 403  km /sec, 

S  wave  velocity:  3  =  2.5A2  km /sec, 

■5 

Density  o  =  2.661  gm/cm"'. 

This  material  is  intended  to  represent  NTS  granite.  A  similar 
series  of  calculations  is  described  by  Perl,  et  a]_.  (1979)  and  Perl 
and  Trulio  (1979).  The  computational  method  and  constitutive  model 
for  the  NTS  granite  is  described  in  these  reoorts. 

The  seven  calculations  were  done  in  a  cylindrical 
(axisymmetric)  geometry  and  differ  only  in  the  source  depth.  These 
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ATI 

Identifier 

Oepth  (m) 

Number  of 
Monitoring 
Stations 

4701 

159.4 

64 

4702 

207.2 

66 

4703 

253.0 

66 

5127 

398.3 

66 

5128 

531.3 

64 

5129 

797.0 

64 

5130 

1000.0 

63 

All  calculations  were  carried  out  to  2.5  seconds. 

The  source  data  orovided  to  S-Cubed  were  the  time  histories  of 
the  "excess"  tractions  ( i . e . ,  tractions  relative  to  their  initial 
equilibrium  values),  7^  (r,z,t),  and  velocities, 

uz>  Cr  (r,z,t),  at  stations  on  a  surface  surrounding  the 
region  of  inelastic  material  response.  The  number  of  monitoring 
stations  for  each  calculation  is  listed  above.  The  position  of 
these  stations  for  a  typical  calculation  is  shown  in  Figure  1. 

Also  provided  to  S-Cubed  were  the  area  of  the  surface  jpon 
which  the  tractions  were  assumed  to  act  and  the  components  of  the 
unit  vector  (nr,  normal  to  that  area.  All  of  these 

quantities  are  necessary  for  calculating  seismograms  in  the  far- 
field,  as  described  in  Section  III  and  Appendix  A. 

In  Figure  2  we  show  the  time  histories  of  the  monitored 
quantities  at  several  oositions  on  the  monitoring  surface  for  ATI 
Number  5127.  The  motions  are  not  entirely  stopped  at  2.5  seconds 
when  the  calculations  were  terminated,  -'owever,  they  are  small  at 
this  time  and  the  decision  to  terminate  the  calculation  seems 
reasonable. 

In  the  Appendix  A,  Section  A. 5,  we  discuss  the  behavior  of  the 
surface  wave  solution  for  periods  that  are  large  compared  to 
duration  of  the  source  calculation.  The  most  important  data  for 
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Figure  2.  The  monitored  velocities  and  tractions  are  plotted  at  selected 
positions  (given  in  terms  of  their  coordinates  in  kilometers) 
for  ATI  s5127.  As  will  be  explained  in  Section  2.4,  the  AH 
calculations  have  positive  z  up  and  minus  signs  in  the  stress- 
strain  relation.  To  be  consistent  with  the  theory  of  Appendix 
A,  it  is  necessary  to  change  the  signs  of  the  vertical  veloc¬ 
ity  and  radial  traction.  We  also  show  the  values  of  the  trac¬ 
tions  at  an  imposed  final  time  steo. 
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Figure  2.  (Continued) 
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this  solution  are  the  static  values  of  TV  ,  T  and  u' 
These  static  tractions  appear  to  be  small  compared  to  the  peak 
values  and  are  probably  difficult  to  compute  accurately  for  this 
reason. 


In  Figure  3  we  plot  the  values  of  the  monitored  solutions  at 
the  last  time  point  versus  position.  Since  the  velocities  should  be 
nearly  zero  at  this  time,  their  angular  dependence  may  not  be  too 
significant.  However,  the  angular  variation  of  the  tractions  is 
important.  Their  raoid  oscillation  with  position  suggests  that 
numerical  errors  may  be  present  at  these  low  stress  levels.  This 
cou-ld  lead  to  significant  errors  in  the  long  period  solutions,  so  we 
will  be  checking  this  part  of  the  solution  with  care. 


2.3  THE  VERTICAL  FORCE  AND  IMPULSE 

An  important  indicator  of  the  numerical  accuracy  of  the 
calculations  is  the  total  vertical  force  and  impulse: 


Fz(t) 


JV 


(r,z,t)  dA 


I 


z 


Fz(')  dt. 


(1) 


where  S  is  the  monitoring  surface.  The  four  deepest  explosions  were 
fully  contained  (n.  3er1,  personal  communication).  Tnen  for  these 
calculations  there  are  no  external  *orces  or  body  forces  influencing 
the  monitored  solutions.  In  this  case,  conservation  of  linear  mo¬ 
mentum  reauires  that  the  force  and  imoulse  vanish  at  late  time. 
Each  of  the  shallowest  three  calculations  formed  a  crater.  But  the 
ejected  material  must  return  to  the  surface  at  late  time  and,  again, 
linear  momentum  will  be  conserved,  though  the  calculation  miqnt  not 
be  run  long  enough  to  satisfy  this  condition. 

In  figure  a  we  show  the  vertical  force  and  imoulse  for  the 
seven  ATI  calculations.  In  each  case  these  Quantities  are  non-zero 
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Figure  3. 


The  monitored  velocities  and  tractions  for  ATI  #5127  are 
olotted  versus  position  for  the  last  time  cycle.  The  moni¬ 
toring  station  positions  are  shown  in  Figure  1  and  are  olotted 
from  r  »  0,  z  =  2.1  at  the  left  to  r  *  2.1,  z  =  0  at  the  riqht. 
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Figure  4.  (Continued) 
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at  the  last  time  step  (at  2.5  seconds),  though  they  are  small,  ■'he 
shallow  calculations  that  cratered  show  no  systematic  differences 
from  the  deeper  contained  source  calculations. 

In  Section  A. 5  of  Appendix  A  we  discuss  the  low  frequency 
asymptotic  behavior  of  the  surface  wave  solution.  If  F  (t)  and 
I2(t)  are  finite,  the  dominant  term  in  the  solution  depends  on  the 
static  value  of  the  vertical  force.  If  the  force  aporoaches  zero  in 
such  a  way  that  the  impulse  remains  non-zero,  the  vertical  force 
term  is  of  the  same  order  as  other  terms  in  the  low  frequency 
asymptotic  solution  and  so  will  make  an  important  contribution. 
However,  if 

lim  F  (t)  =.  lim  I  (t)  -  0  ,  (2) 

t>°°  t>ao 

the  vertical  force  term  is  a  higher  order  term  that  makes  little 
contribution  at  long  periods.  Conservation  of  linear  momentum 
requires  that  (2)  be  satisfied,  and  our  asymptotic  analysis  shows 
that  this  is  an  important  constraint  on  the  long  period  solution. 

2.4  CORRECTING  THE  DATA  TO  ZERO  THE  VERTICAL  FORCE  AND  IMPULSE 

M 

A  straightforward  procedure  may  be  devised  to  correct  the  T 
to  satisfy  the  conditions  (2).  Basically,  we  assume  that  the  final 
values  of  the  force  and  impulse  are  nonzero,  either  because  there 
are  some  numerical  errors,  or  because  the  calculation  is  not  auite 
finished  (i.e.,  some  small  motions  are  still  propagating  through  the 
monitoring  surface).  We  then  add  an  additional  time  step  to  correct 
the  monitored  solutions  to  satisfy  (2).  Tnis  is  an  ad  hoc  procedure 
that  is  least  troublesome  if  the  correction  is  "small";  that  is,  if 
we  start  with  calculations  where  (2)  is  nearly  satisfied,  as  is  the 
case  here. 

Our  procedure  for  "correcting"  the  monitored  tractions  is  as 
follows.  Let  t  «  T^  at  the  final  time  step  of  the  calculation. 
Then  assume  that 
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In  practice,  we  will  always  have  K  *  0  and  L  ±  0. 
another  time  step,  at  =  such  that 

Fz  (T2)  *  I2  (T2)  =  0. 


Then  we  take 


This  is  done  by  setting 


T2  <V  -  T2  <Tl)  *  sTZ 
at  each  point  on  the  monitoring  surface. 


From  (1)  we  see  that  the  vertical  force  will  be  zero  if  we  choose 
51^  to  satisfy 


sTz  dA  .  -  K  . 


If  K  and  i  have  the  same  sign,  we  can  zero  the  impulse  by  an 
appropriate  choice  for  at.  Let  Fz(t)  be  linear  between  and 
T2.  That  is, 

Fz(t)  -  K  -  (t  -  Tj)  K/at  (7) 

for  ti.  [Tp  T^].  Then 

at  .  (3) 


leads  to 


!r  '<  and  L  have  the  oooosite  sign,  a  more 


complex  functional  form  could  be  chosen  cor  "2(t)  between  t  and 
Tp  to  zero  the  impulse.  However,  in  most  cases  we  prefer  to 
retreat  to  an  earlier  time  step  where  K  and  l  have  the  same  sign  and 
then  aoply  the  procedure  outlined. 

We  have  not  yet  shown  how  to  select  the  af  to  satisfy  (F) 
while  minimizing  the  perturbation  to  the  individual  7  time 
histories.  Tnere  are  other  possibilities,  but  we  have  chosen  to 
minimize 
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while  requiring  that  (6) 
minimization  problem  which 
multiplier.  The  solution  is 


be  satisfied.  This  is  a 
can  be  solved  using  a 


classic 
Lagrange 


In  summary,  the  tractions  and  displacements  are  comDuted  to 
time  T,.  Then  another  time  step,  at  from  (8),  is  added  to  take 
the  calculation  to  a  new  final  time  T».  The  vertical  tractions  at 

^  M 

T^  are  computed  from  (5),  using  (10)  to  compute  5rz.  The 
other  quantities  are  extended  as  follows: 

<V  (T2)  =  uj!  (Tj)  , 

uz  *  uz  <T1>  *  (12) 

rr  (T2)  »  T?  (V  * 

In  Table  1  we  list  the  parameters  used  for  correcting  the 
vertical  force  and  impulse  from  the  ATI  calculations.  Mote  that  fo r 
two  of  the  calculations  (4702  and  5130)  some  of  the  data  provided 
were  discarded  and  the  processed  data  were  terminated  at  a  time  when 
the  force  and  impulse  had  opposite  signs. 

The  important  numbers  in  Table  1  are  KIF  and  at.  We  see  that 
the  imposed  final  time  step  is  not  too  large.  Calculation  5127  is 
typical  and  the  correction  is  indicated  in  Figure  5.  The  K/T  give 
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TA8LE  1 


PARAMETERS  FOR  CORRECTING  THE  VERTICAL  FORCE  AND  IMPULSE  TO  ZERO 


ATI 

Identifier 

Depth 

T1 

(sec) 

<VTi>  MTi> 

(1017  (1016 
dynes) (dyne-sec) 

r» 

do17 

dynes) 

k  /: 

At 

(sec) 

4701 

159.4 

2.50 

4.5 

-7.8 

7.2 

0.62 

0.35 

4702* 

207.2 

2.45 

2.6 

-1.5 

5.4 

0.49 

0.12 

4703 

253.0 

2.50 

-1.9 

2.9 

7.6 

-0.25 

0.30 

5127 

398.5 

2.50 

-4.1 

8.3 

8.1 

-0.50 

0.41 

5128* 

531.3 

2.37 

-2.2 

2.8 

8.5 

-0.26 

0.25 

5129 

797.0 

2.50 

-1.9 

3.4 

4.5 

-0.42 

0.36 

5130 

1000.0 

2.50 

-0.46 

14.4 

5.8 

-0.08 

0.63 

*  The  last  time  step  used  for  subsequent  calculations  is  earlier 
than  the  last  time  step  provided. 
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the  ratio  of  the  traction  correction  (aT")  to  the  computed  value 
at  the  last  time  step,  see  (10).  For  all  but  number  5120,  this 
correction  is  rather  large.  For  example,  for  Number  5127  we  have 


0.50  Tj 


That  is,  to  zero  the  vertical  force  and  impulse,  we  impose  an 
“upward"  traction  at  each  node  that  is  half  the  size  of  the  computed 
final  value.  The  values  of  the  tractions  used  in  the  calculations 
(including  the  imposed  time  step)  were  plotted  in  Figure  2.  In  many 
cases,  the  corrections  are  small  in  absolute  terms  because  the  final 
computed  values  were  small  to  start  with. 

How  important  are  these  corrections?  First,  they  mainly 
affect  long  period  seismic  waves.  That  is,  the  surface  waves.  Our 
asymptotic  analysis  in  Appendix  A  shows  that  the  contribution  of  the 
vertical  traction  to  long  period  surface  waves  is  rather  small  as 
long  as  the  force  and  imoulse  do  vanish.  For  these  reasons,  we  do 
not  expect  the  details  of  the  correction  to  be  very  important.  We 
will  later  show  in  Section  3.4  that  this  is  indeed  the  case. 


SYSTEMS,  SCIENCE  AND  SOFTWARE 


! 


III.  Ms  and  mb  Estimates  for  ATI  Granite  Calculations 
3.1  INTRODUCTION 

The  ATI  source  calculations  described  in  the  previous  section 
were  processed  with  the  methods  described  in  Appendix  A  to  compute 
synthetic  body  and  surface  wave  seismograms  in  realistic  earth 
models.  From  these  seismograms  conventional  m^  and  M$  were 
determined.  In  this  section  we  describe  these  synthetic  seismogram 
calculations  and  analyze  the  results. 

The  basic  equation  for  the  seismogram  synthesis  is  (Section 

A. 2). 


*  -/H 


'*t^ 

•'  j 


-i'  M 
-  S  *u  . 
J 


dA 


(14) 


i  i  m 

where  G-  and  SM,  are  Green's  functions  while  r  and 

u .  are  the  monitored  tractions  and  disdacements.  For  most  of 

the  development  in  Appendix  A,  the  tractions  are  replaced  by  the 

actual  stress  components.  In  that  form  the  basic  equation  for 

axisymmetric  sources  is  (4.5).  In  terms  of  tractions,  and  a  general 

monitoring  surface  like  that  in  Figure  1,  the  vertical  disolacements 

are 


-/ 


Gr  fr 


*s"  Tz 


dA  + 


/[ 


M  ,  -Z  M 
S  u  n  +  S  u  n 
rr  r  r  zr  z  r 


dA 


(15) 


The  G1.  and  sL  are  the  azimuthally  averaged  Green's 

J  mm 

functions  given  in  Apoendix  A.  The  T. ,  u.  and  n,  are 

J  J  J 

the  quantities  provided  by  ATI.  The  quadratures  reduce  to  a  sum 
over  the  nodes  with  the  dA  being  the  aoprooriate  incremental  area. 
Equation  (15)  is  used  for  all  calculations  in  this  section.  The 
difference  between  body  and  surface  waves  is  in  the  Green's 
functions. 
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3.2  FAR-FIELD  RAYLEIGH  WAVES  !M  A  HAl.FSPACE 

We  begin  our  study  of  the  Rayleigh  waves  excited  by  the  atj 
granite  calculations  by  examining  the  Rayleigh  oulse  in  a  half¬ 
space.  The  waveform  is  simple  and  unusual  features  are  easy  to 
see.  In  dispersive  real  earth  models,  peculiarities  due  to  the 
source  are  relatively  difficult  to  identify. 

The  halfspace  is  that  of  the  source  calculations  fa  *  4,ao3 

•j 

'<m/sec,  8  =  2.542  km/sec,  o  =  2.661  gm/cmJ).  The  halfsoace 
Rayleigh  pulses  were  computed  at  1000  kilometers  and  were  filtered 
by  a  WWSSN  15-100  long  oeriod  seismometer  resoonse.  The  effects  of 
attenuation  (0  *  10  in  the  top  kilometer  and  300  elsewhere)  were 
included  by  multiplying  the  Rayleigh  wave  soectrum  by  the  following 
factors: 


Period 

Attenuate 

(sec) 

Factor 

25.0 

0.73 

20.0 

0.65 

14.0 

0.48 

10.0 

0.31 

6.25 

0.097 

5.0 

0.0006 

The  halfspace  Rayleigh  oulse  is  shown  in  Figure  6  compared  to 
the  Rayleigh  oulse  from  an  ROP  source  at  the  same  depth  as  the  att 
explosion.  The  ROP  is  that  of  Mueller  and  Murphy  (1971)  for  granite 
at  a  depth  of  398.5  meters,  ’here  is  little  perceptible  difference 
between  the  Rayleigh  waves  from  the  ATI  sources  ano  the  R0-  source. 
We  do  see  a  little  long  period  "noise1*  on  several  of  the  seismograms 
for  the  ATI  sources,  but  generally  they  look  quite  good. 
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?eak-to-Peak 

Source  Amplitude  Depth 

Identifier  (Microns)  (m) 


Peak-to-Peak 

Source  Amplitude  Depth 

Identifier  (Microns)  (m) 


3.3  CONTRIBUTION  OF  VERTICAL  TRACTIONS 

As  was  discussed  in  Section  2. a,  an  ad  hoc  correction  was 
introduced  into  the  vertical  tractions  to  force  them  to  satisfy  the 
required  condition  that  the  total  downward  force  and  impulse  vanish 
at  late  time.  The  influence  of  this  correction  on  the  contribution 
of  the  vertical  traction  terms  to  the  solution  should  become  larger 
with  increasing  period.  We  do  not  expect  ‘the  correction  to  have 
much  influence  on  periods  below  two  seconds  that  are  important  for 
mb<  However,  the  important  periods  for  surface  waves  are  much 
larger  than  two  seconds. 

In  Figure  7  we  compare  the  halfspace  Rayleigh  wave  solution 
from  Figure  6  to  the  Rayleigh  wave  computed  with  only  the  vertical 
tractions.  The  latter  were  computed  with  an  incorrect  sign,  so  the 
actual  vertical  traction  Rayleigh  wave  is  inverted  compared  to  that 
in  the  figure.  That  is,  it  has  the  same  polarity  as  the  total 
solution.  The  Rayleigh  wave  spectra  for  these  seismograms  are 
plotted  in  Figure  8. 

The  dominant  period  of  the  total  halfspace  Rayleigh  wave  is 
about  nine  seconds.  cor  the  vertical  traction  component  the 

dominant  period  is  slightly  shorter,  about  seven  seconds.  rne 
peak-to-peak  amplitudes  and  their  ratios  shown  in  figure  7  are  not 
corrected  for  the  small  difference  in  instrument  response  at  these 
periods.  If  included,  this  correction  would  increase  the  ratios  by 
about  15%. 

For  M  our  interest  is  in  periods  longer  than  nine  seconds. 

The  soectral  comparison  in  Figure  3  shows  th'*  the  vertical  traction 

contrioution  is  less  for  longer  oeriods  than  it  is  at  the  period 

dominating  the  halfsoace  Rayleigh  waves. 

The  ratios  in  Figure  7  and  the  spectral  comparisons  indicate 
that  the  contribution  of  the  vertical  tractions  decreases  with 
source  depth.  For  the  three  shallowest  depths,  the  vertical  trac¬ 
tion  contribution  is  within  a  factor  of  two  or  three  of  the  total 
solution  for  M  periods.  This  is  merely  a  suggestion  that  we 
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Figure  7.  Comparison  of  halfspace  Rayleigh  waves  for  the  total  solution 
to  those  computed  from  only  the  vertical  tractions.  The 
latter  are  inverted.  The  peak-to-peak  amplitudes  are  shown 
for  each  seismogram  and  the  ratio  is  given  for  eacri  pair. 
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FREQUENCY  (HZ) 


Figure  3 


The  Rayleigh  wave  spectra  for  the  seismograms  of  Figure  7 
are  plotted.  These  spectra  do  not  include  the  effect  of 
0  or  the  seismometer.  The  smoother  spectrum  that  is 
larger  over  most  of  the  freauency  range  is  the  total 
spectrum. 
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should  have  less  confidence  in  the  solution  at  the  shallower 
depths.  Of  course,  we  have  no  way  to  estimate  the  relative 

importance  of  errors  from  the  ad  hoc  correction  to  the  vertical 

tractions  (which,  after  all,  forces  them  to  satisfy  a  global 
conservation  law),  compared  to  other  errors  in  the  calculations. 

3.4  AN  EQUIVALENT  RDP  FRCM  SURFACE  WAVES 

Toe  standard  way  to  represent  the  source  for  underground 
nuclear  exDlosions  is  in  terms  of  a  spherically  symmetric  point 
source.  An  important  motivation  for  doing  multi-dimensional  source 
calculations  is  to  study  the  influence  of  effects  that  cannot  be 
represented  by  a  one-dimensional  source.  One  way  to  do  this  is  to 
compare  seismograms  computed  from  the  complex  source  with  those  from 
a  one-dimensional  (reduced  displacement  ootential  or  ROP)  source. 

Such  a  comparison  was  made  in  the  time  domain  in  Figure  6. 

In  Appendix  3  we  point  out  that  an  "eouivalent  RDP"  repre¬ 
sentation  can  be  computed  for  the  two-dimensional  source.  This 

equivalent  source,  called  '?  ,  can  be  comouted  from  the  equation 


where  WA(,  is  the  spectrum  of  the  Rayleigh  waves  from  the 
two-dimensional  sources  in  Figure  6.  The  is  the  spectrum  of 
the  corresponding  seismogram  computed  with  the  reduced  velocity 
ootential,  ?  . 

Another  way  to  look  at  the  'ia  is  that  it  is  the  source  that 
we  would  deduce  from  Rayleigh  wave  recordings  of  an  explosion,  if  we 
had  perfect  knowledge  of  the  path  and  if  we  assumeo  that  the  source 
could  be  represented  by  an  RDP. 

The  semi-empirical  Muel ler/^urphy  model  (Mueller  and  Murphy, 
1971;  Murphy,  1977)  predicts  that  the  long  oeriod  source  level  (’Fw) 
depends  on  depth  according  to 
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for  explosions  at  a  fixed  yield.  The  Mueller /Murphy  granite  source 
functions  for  the  seven  depths  of  the  150  KT  ATI  calculations  are 
shown  in  Figure  9. 

The  "equivalent  RDP"  reoresentations  for  the  ATI  calculations 
were  computed  with  (16)  using  the  spectra  of  the  seismograms  in 
Figure  6.  These  are  shown  in  Figure  10.  The  amplitude  scales  vary 
from  case  to  case,  so  a  Mueller /Murphy  source  function  is  plotted 
with  each  case  for  reference. 

/}. 

The  high  frequency  behavior  of  the  in  Figure  10  is  not 
what  we  expect  for  a'true  RDP  source.  We  see  that  the  amplitudes  go 
off-scale.  Mechanically,  this  is  because  the  Rayleigh  wave  spectra 
for  the  ROP  sources,  Wp  ,  fall  off  faster  than  u>“^  at  high 
frequencies,  while  the  spectra  for  the  two-dimensional  sources, 
WAC  (Figure  8),  are  not  even  monotonically  decreasing  at  high 
frequency.  Part  of  this  may  be  due  to  numerical  error,  since  the 
surface  wave  synthesis  method  may  be  losing  accuracy  at  the  high 
frequencies.  However,  there  are  also  plausible  physical  expla¬ 
nations  for  the  observed  behavior  of  the  WAC  and  1^1  .  First, 
the  free  surface  reflections  in  the  inelastic  ATI  calculations  are 
not  perfect  scaled  replicas  of  the  direct  waves,  as  they  are  in  the 
calculations  of  W^.  Second,  the  ATT  calculations  presumably 
include  some  contribution  from  spall  closure.  These  physical  expla¬ 
nations  for  the  high  frequency  character  of  the  1’?  I  are  supoorted 
by  the  way  the  soectral  holes  shift  to  lower  freauency  with  in¬ 
creasing  depth. 

3.5  PAR-FIELD  RAYLEIGH  WAVES  VI  A  REALISTIC  cART’-i  MODEL 

The  fundamental  mode  Rayleigh  waves  for  the  ATI  source 
calculations  were  computed  for  a  crust  and  upper  mantle  model  for 
the  central  United  States.  The  velocity  model  is  essentially  that 
given  by  McEv illy  (196A)  while  the  Q  model  is  roughly  compatible 
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FREQUENCY  (Hz) 


Figure  10.  The  |fei  from  the  seismograms  in  Figure  6 
The  amplitude  axes  are  in  104  m3. 

Murphy  I  f  I  for  a  depth  of  700  meters  (see 
shown  for  reference. 
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FREQUENCY  (Hz) 


iqure  10.  (continued) 
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with  the  Rayleigh  wave  attenuation  factors  of  Trygovason  (1965). 
These  models  are  listed  in  Table  2.  To  account  for  the  mismatch 
between  the  near-surface  prooerties  in  the  source  region  and  those 
in  the  propagation  model,  we  used  the  techniaue  of  Bache,  Rodi  and 
Harkrider  (1978),  which  is  described  in  Section  A. 4.  With  this 
method  separate  models  are  used  for  the  source  region  and 
propagation  path.  The  two  differ  only  in  the  top  three  kilometers. 

The  computed  Rayleigh  waves  are  shown  in  Figure  11.  Trie 
structure  has  a  strong  Airy  phase  near  16  seconds.  The  somewhat 
strange  appearance  of  the  waveform  results  from  the  simultaneous 
arival  of  5  to  3  second  waves  from  the  inverse  branch  with  the  20  to 
22  second  waves. 

The  amplitude  for  computing  Mg  was  taken  from  the  largest 
trough  to  the  following  peak.  The  period  of  this  cycle  is  about  14 
seconds.  The  M$  formula  of  Marshall  and  Basham  (1972)  was  used. 
At  3000  kilometers  this  is 

Ms  »  log  A  +  1.38  +  P(T)  (18) 

where  A  is  the  maximum  zero-to-peak  amplitude  in  millimicrons,  P(T) 
is  a  period-dependent  path  correction  and  1.38  is  the  distance 
correction. 

The  M$  values  are  listed  in  Table  3  and  are  plotted  versus 
source  depth  in  Figure  12.  As  another  display  of  the  depth-depen¬ 
dence  of  the  Rayleigh  wave  excitation,  we  also  plot  the  logarithms 
of  the  Rayleigh  wave  spectral  amplitudes  at  20,  14  and  10  seconds. 
The  depth  deoendence  is  greatest  at  20  seconds  and  decreases  at 
shorter  periods.  The  M$  depth  deoendence  lies  between  those  for 
the  14  and  20  second  spectral  amplitudes. 

A  simpler  display  that  is  independent  of  the  propagation  model 
is  to  plot  the  logarithm  of  the  spectral  amplitude  for  the  halfspace 
Rayleigh  waves  shown  in  Figure  7.  These  are  shown  in  the  figure  for 
20,  10,  5  and  1  second  periods.  The  20  and  10  second  curves  are 
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CRUSTAL  MODELS  FOR  SUReAC£  WAVE  CALCULATIONS 


Oeoth 

Thickness 

a 

3 

a 

(km) 

(km) 

(km/sec) 

(km/sec 

(gm/crT) 

Source 

Region  Mooel 

1 

1 

4.403 

2.542 

2.661 

2 

1 

4.403 

2.542 

2.661 

3 

1 

4.403 

2.542 

2.551 

11 

8 

6.1 

3.5 

2.7 

20 

9 

6.4 

3.58 

2.9 

38 

18 

6.7 

3.94 

2.9 

62 

24 

8.15 

4.75 

3.3 

102 

AO 

8.2 

4.51 

3.3 

120 

18 

8.7 

4.80 

3.6 

8.7 

4.80 

3.6 

Propagation  °ath  Model 

Extend  Layer  a  (o  =6.1  km/sec)  to  the  surface 


Figure  11.  Rayleigh  waves  in  a  shield  crustal  model  for  seven  ATI  granite 
calculations.  The  range  is  3000  km  and  the  WWSSN  15  to  100 
instrument  response  was  included.  Zero  on  the  time  scale  cor¬ 
responds  to  750  seconds  after  source  initiation. 
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TABLE  3 


Calculation 


Ms  FOR  SEVEN  ATI  CALCULATIONS 
OF  EXPLOSIONS  IN  GRANITE 


Depth  A 

(km)  (microns)  ( 


nearly  the  same  as  those  for  the  disoersive  crustal  model,  as 
expected.  We  see  that  the  dependence  on  depth  continues  to  change 
as  we  go  to  shorter  periods. 


The  spectral  values  in  Figure  12  enhance  our  confidence  that 
the  M$  computed  with  the  particular  crustal  model  chosen  is  a  true 
representation  of  the  surface  waves  from  the  source  calculations. 
The  increases  slightly  as  the  depth  goes  from  159  to  207 
meters.  There  is  a  sharp  decrease  at  253  meters  followed  by  an 
increase  as  depth  increases  to  531  meters.  The  Mg  values  then 
fall  off  for  increasing  depth. 

Another  interesting  comparison  is  shown  in  Figure  13.  In  this 
plot,  we  show  the  Mg  and  20  second  spectral  amplitude  (WA(0  from 
Fiqure  12  compared  to  the  20  second  value  of  the  equivalent  RD? 

$  (  tg)  from  Figure  10.  The  trend  for  all  three  measures  of  the 

surface  wave  excitation  is  the  same. 

Also  shown  in  Figure  13  is  the  20  second  spectral  amplitude  of 
the  Muel ler /Murphy  P.OP  from  Figure  9.  This  is  nearly  the  expected 
|  depth  dependence  of  M$  predicted  by  the  Mueller /Murphy  model  (the 

small  dependence  of  the  Rayleigh  wave  excitation  on  depth  is 
ignored).  At  greater  depths,  where  the  interaction  with  the  free 
surface  may  be  nearly  elastic,  the  depth  dependence  of  the  Rayleigh 
wave  excitation  for  the  ATI  calculations  is  not  too  different  from 
that  for  the  Mueller /Murphy  model.  However,  at  shallow  depths  the 
two-dimensional  calculations  predict  a  much  different  depth 
dependence. 

3.6  500Y  WAVE  MAGNITUDE 

Synthetic  body  wave  seismograms  were  computed  using  techniaues 
similar  to  those  for  computing  the  Rayleigh  waves.  The  basic 
■  equation  is  (15)  with  Green's  functions  appropriate  for  teleseismic 

body  waves. 

The  body  wave  seismograms  were  for  a  range  of  A000  km.  The 
earth  model  included  detailed  crustal  structures  for  the  source  and 
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receiver  regions  (Table  4).  Rather  than  use  a  detailed  upper  mantle 

model,  we  represented  this  part  of  the  path  by  a  constant  geometric 

spreading  factor,  which  is  reasonable  for  this  range.  For  all  the 

seismograms  shown  in  this  section  (as  well  as  for  the  S-Cubed 

5  1 

sources  discussed  in  Section  V),  we  used  5.6  x  10  km  for  the 

effective  1/R.  This  value  is  somewhat  arbitrary.  Langston  and 

Helmberger  (1975)  give  about  7.6  x  10“^  km-*  for  this  range  in  a 

Jeffreys-Bullen  earth  model.  We  chose  the  ray  parameter  (source 

takeoff  angle)  to  be  appropriate  for  a  range  of  4000  km  in  the  model 

HWNE  (Helmberger  and  Wiggins,  1971).  Later  calculations  indicated 
5  1 

that  6.4  x  10  km  is  the  best  estimate  for  the  effective  1/R 
at  this  range  in  HNME.  The  seismograms  are  easily  scaled  for  a  1/R 
different  than  the  value  we  used. 

The  ray  parameter  for  the  body  wave  synthesis  corresponds  to  a 
source  takeoff  angle  of  20.3  degrees.  The  calculations  are  not  very 
sensitive  to  tnis  choice.  It  is  also  necessary  to  account  for  Q. 
This  was  done  by  multiplying  the  spectrum  by  the  causal  Q  operator 
(Strick,  1970)' 


where  f  is  frequency  in  Hertz.  Our  calculations  were  done  with  t*  * 
0.3.  This  is  an  important  path  dependent  parameter  which  is  much 
discussed  and  debated  in  the  literature.  We  believe  appropriate 
values  for  short  period  0  wave  synthesis  range  from  0.5  to  1.2,  so 
our  choice  lies  between.  Finally,  the  seismograms  were  filtered  by 
the  response  of  the  KS36000  seismometer.  The  amplitude  response  is 
plotted  in  figure  14. 

The  synthetic  body  wave  seismograms  are  plotted  in  Figure  15. 

Two  amplitude  measurements  were  made  from  each  record,  one  on  the 

"b"  phase  (first  peak  to  first  trough)  and  one  on  the  "c"  phase 

b 

(first  trough  to  second  peak).  Body  wave  magnitudes,  m^  and 
m5,  were  computed  from  these  amplitudes  using 
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TABLE  4 

CRUSTAL  MODELS  FOR  BODY  WAVE  CALCULATIONS 


Depth 

(km) 

Thickness 

(km) 

a 

(km/sec) 

S 

(km/ sec) 

(gm/cm3) 

3.0 

3.0 

SOURCE  REGION 

4.403 

2.542 

2.661 

8.0 

5.0 

5.35 

2.79 

2.70 

20.0 

12.00 

6.0 

3.5 

2.70 

1.70 

1.70 

RECEIVER  REGION 

4.0 

2.31 

2.3 

3.00 

1.30 

5.1 

2.94 

2.5 

20.00 

17.00 

6.0 

3.5 

2.8 
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Figure  15,  Synthetic  body  wave  seismograms  for  seven  ATI  calculations  of 
explosions  in  granite.  The  maximum  oeak-to-peak  amplitude  m 
microns  at  one  Hertz  ’s  listed  at  the  left  of  eacn  seismogram. 
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m£  =  log  |y]+  3.25,  U?) 

where  A  and  T  are  the  appropriate  amplitude  (corrected  for 
instrument  response)  and  period  and  3.25  is  the  distance  correction 
for  4000  kilometers. 

The  magnitudes  are  listed  in  Table  5.  Mote  that  the  apparent 
period  of  the  b  and  c  phases  decreases  slightly  with  depth  to  253 
meters  (the  maximum  depth  at  which  cratering  occurred),  then 
increases  thereafter. 

The  m£  and  are  plotted  versus  depth  in  Figure 

16.  Also  plotted  are  the  M  from  Figure  8  and  the  residuals, 

h  r  ^ 

mP  -  and  mP  -  M„.  The  dependence  of  m.  and  M„ 

b  s  b  s  d  s 

on  depth  is  auite  similar  except  for  the  point  at  253  meters  depth. 
At  that  point  there  is  a  sharp  minimum  in  the  M  that  is  not 
present  in  the  body  wave  magnitudes. 


3.7  ANALYSIS  OF  mfa  DATA 

The  short  period  P  wave  seismograms  for  the  two-dimensional 
ATI  calculations  can  be  studied  by  comparing  them  to  analogous 
seismograms  computed  with  a  one-dimensional  RDP  source.  Such  a 
comparison  is  shown  in  Figure  17.  The  seismograms  for  the 
two-dimensional  ATI  sources  from  Figure  15  appear  in  the  center 
column.  In  the  left  column  we  show  seismograms  at  the  same  depth 
for  a  constant  RDP  source  which  is  the  Mueller /Murphy  RDP  computed 
at  the  shallowest  depth,  159.4  meters.  Thus,  only  the  P-pP  lag  time 
changes  with  deoth  for  the  seismograms  in  the  left  column. 

In  the  right  column  in  Figure  17  we  show  the  seismograms  for 
the  Muel ler /Murphy  RDP  source  with  the  depth  scaling  included.  The 
source  functions  used  for  these  seismograms  were  plotted  in  Figure 
9.  With  increasing  depth,  the  source  function  gets  smaller  and  the 
peak  moves  to  higher  freouencies.  The  b  phase  magnitude,  m^, 
and  the  associated  period,  Tb,  are  listed  with  each  seismogram  and 
they  have  the  expected  trend. 
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a.  ?CR  SZVZN  ATI  CALCULATIONS  C?  INRLCSICNS  IN  GRANITE 
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3127 

398.5 

134 

0.73 

5.52 

244 

0.34 

5.71 

5.23 

531.3 

245 

0.36 

5.71 

347 

0.39 

3.34 

5129 

797.0 

219 

0.92 

5 . 62 

340 

0.97 

5.79 

5130 

1000.0 

202 

0.93 

5.59 

325 

1.01 

5.75 

CONSTANT  R DP  TWO  -  D  E1ENS I ONAL  MCELLER/MUSPHY  RD? 

SOURCE 

Figure  17.  Comparison  of  the  short  period  body  wave  seismogram  for  the 
ATI  granite  source  calculations  with  seismograms  computed  with 
a  spherically  symmetric  RDP  source. 
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Comparison  of  the  waveforms  in  Figure  17  suggests  that  the 
two-dimensional  source  calculations  give  body  waves  that  do  look 
much  like  the  simpler  P  +  pP  seismograms  from  the  one-dimensional 
sources.  We  see  that  the  second  peak  is  initially  the  same  size  as 
the  first,  then  gets  larger  as  the  P  and  pP  constructively  interfere. 

The  ATI  calculations  include  nonlinear  material  behavior  all 
the  way  to  the  free  surface.  In  fact,  the  first  three  formed  a 
crater.  Thus,  we  would  expect  the  free  surface  reflected  phase  to 
be  somewhat  different  from  the  pP  of  elastic  theory. 

In  Figure  18  we  again  plot  the  body  wave  seismograms  for  the 
ATI  calculations  and  compare  them  to  seismograms  computed  with  the 
Mueller /Murphy  RDP.  This  time  the  source  is  modified  to  completely 
suppress  the  pP  (left  column)  or  to  halve  it.  The  waveform 
comparison  suggests  that  substantial  suppression  of  pP  is  not 
occurring  for  the  greater  source  depths.  However,  this  may  be  a 
more  reasonable  model  at  the  shallow  depths. 

These  seismogram  comparisons  are  placed  on  a  more  quantitative 
basis  in  Figures  19,  20  and  21.  In  these  figures,  we  compare  the 
magnitude  and  period  data  taken  from  the  seismograms  in  Figures  17 
and  18. 

In  Figures  19  and  20  we  plot  the  magnitude  data  from  the 
seismograms.  First,  we  point  out  that  we  are  most  concerned  with 
comparing  the  trends  with  increasing  depth.  However,  the  absolute 
magnitudes  are  quite  close,  especially  for  depths  of  350  to  550 
meters.  This  is  interesting,  but  we  are  not  quite  sure  what  it 
means. 

From  the  comparisons  in  Figures  19  through  21,  we  draw  the 
following  conclusions: 

•  When  the  RDP  source  calculations  include  an 
elastic  pP,  as  in  Figure  17,  the  magnitude 
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Fiaure  18.  The  body  wave  seismograms  Tor  the  ATI  sources  are  comoared  to 
'  ”  ""  those  for  modified  versions  of  the  Mueller/  Murphy  RDP. 
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l*  plotted  versus  source  depth. 
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increases  with  depth.  The  ATI  2-0  calculations 
do  not  show  this  behavior,  especially  at 
shallow  depths. 

•  The  closest  agreement  with  the  magnitude  depth 
dependence  of  the  2-D  calculations  is  when  all 
or  part  of  pP  is  suppressed  in  the  RDP  source 
calculations.  The  filter  that  suppresses  pP  is 
certainly  more  complicated  than  the 
frequency-independent  factor  used  here,  but  the 
ATI  calculations  seem  to  support  the  idea  that 
some  suppression  has  taken  place. 

•  Comparison  of  the  measured  periods  shows  the 
worst  agreement  for  the  pP  suppressed  cases  and 
the  best  agreement  for  the  constant  RDP 
source.  This  is  an  indication  that  the  ATI 
calculations  do  not  show  the  large  corner 
frequency  shift  with  depth  predicted  by  the 
Muel ler /Murphy  model. 

3.8  SUMMARY 

We  have  computed  body  and  surface  wave  seismograms  for  the  ATI 

calculations  and  have  determined  the  m.  and  M„  associated  with 

b  s 

these  calculations.  The  actual  values  of  the  magnitudes  are  to  some 
degree  dependent  on  the  details  of  the  path  models  chosen.  However, 
the  changes  with  source  depth  appear  to  be  nearly  path  independent. 

Neither  m^  nor  M  is  very  strongly  dependent  on  depth. 
The  largest  effect  is  on  M  for  the  shallow  depths,  with  the  value 
at  253  meters  depth  being  most  noticeable.  The  source  coupling  into 
both  body  and  surface  waves  decreases  steadily  below  500  meters.  At 
shallower  depths,  nonlinear  free  surface  effects  are  likely  to 
obscure  the  smoother  effect  of  changing  overburden  pressure. 
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The  results  for  the  ATI  calculations  have  been  analyzed  by 
comparing  to  seismograms  computed  with  the  one-dimensional  reduced 
displacement  potential  model  usually  used  to  represent  the  seismic 
source  function.  This  comparison  highlights  the  influence  of 
two-dimensional  effects.  The  standard  for  comparison  is  the 
Muel ler /Murphy  semi -empirical  RDP  for  granite,  which  is  based  upon 
observations  of  ‘SHOAL.  The  magnitudes  for  the  ATI  two-dimensional 
calculations  are  in  close  agreement  with  those  from  the 

Mue  Her  /Murphy  model  at  depths  near  400  meters.  The  long  period 
level  of  the  Mueller /Murphy  model  decreases  more  rapidly  with  depth 
than  does  that  of  the  ATI  calculations,  as  is  seen  in  the  M 
comparison  in  Figure  13.  The  M$  for  the  ATI  calculations  also 
includes  the  effect  of  nonlinear  interaction  with  the  free  surface, 
including  spallation,  which  is  not  easily  quantified. 

The  m^  comparison  between  one-  and  two-dimensional  sources 
was  made  in  the  previous  section  and  our  conclusions  were 
summarized.  Again,  the  depth  dependence  of  both  amplitude  and 
corner  frequency  is  less  than  predicted  by  the  Muel ler /Murphy 
model.  The  comparison  indicates  that  the  pP  is  influenced  by 
nonlinear  interaction  with  the  free  surface. 

In  the  next  two  sections  we  will  be  describing  similar  two- 
dimensional  calculations  done  by  S-Cubed  for  explosions  in  granite. 
Then  in  Section  VII  we  will  discuss  the  entire  set  of 
two-dimensional  granite  calculations  and  summarize  the  important 
conclusions. 


IV.  S— CUBED  SEMITE  CALCULATIONS 


4.1  introduction 

Four  calculations  of  explosions  in  granite  were  done  at 
S-Cubed.  These  were  two-dimensional,  axisymmetric  calculations; 
that  is,  the  geometry  is  essentially  the  same  as  for  the  ATI 
calculations  discussed  in  previous  sections.  In  this  section,  we 
describe  the  general  features  of  these  calculations  and  the  stresses 
and  displacements  monitored  in  the  elastic  response  regime.  These 
data  were  orocessed  with  the  techniques  described  in  Appendix  A  to 
construct  synthetic  seismograms  and  the  results  will  be  described  in 
Section  V. 

The  S-Cubed  calculations  are  intended  to  represent  explosions 
in  the  PILEDRIVER  environment.  The  first  of  these  was  an  attempt  to 
specifically  model  the  PILEDRIVER  event.  Rimer,  et  aj_.  (1979) 
describe  this  calculations  and  compare  computed  velocity  and 
displacement  time  histories  to  the  observed  data  at  some  25 
near-field  gauges.  This  is  essentially  the  entire  near-field  data 
base. 

The  PILEDRIVER  geometry  used  for  the  calculation  and  the  gauge 
locations  are  shown  in  Figure  22.  The  data  were  collected  by  SRI 
(Hoffman  and  Sauer,  1969),  except  the  shot  level  stations  3-SL  and 
16— SL  which  were  collected  by  Sandia  Laboratories  (Perret,  1968). 

The  layered  geologic  model  was  constructed  from  P  wave  arrival  time 
data  in  these  reoorts  and  in  a  description  of  the  nearby  uaRDhat 
event  (Swift,  1962).  The  shear  wave  velocity  was  based  on  an  i 

assumed  °oisson'3  ratio  of  about  0.3. 

As  described  oy  Rimer,  et  aj_.  (  19791  ,  the  S-Cubed  PILEDRIVER  ■'] 

calculation  included  an  effective  stress  law  to  account  for  the  j 

water  present  in  the  jointed  granite  below  the  water  table.  The 

void  crush-up  that  is  important  above  the  water  table  was  accounted  j 

for  by  a  P-a  model  (Cherry,  et  al.,  1975).  Tensile  cracking  is 
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Figure  22.  Assumed  layered  geology  and  location  of  gauge  stations  for  PILEDRIVER  and  other  S-Cubed 


important  in  the  entire  region  surrounding  the  explosion  and 
especially  near  the  surface  where  spallation  and  subseouent  spall 
closure  are  important  contributors  to  the  ground  motions. 


1 

y- 
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In  Figures  23  and  24  we  compare  the  computed  and  observed 
velocity  time  histories  at  two  of  the  stations,  one  at  shot  level 
and  one  on  the  surface.  These,  and  the  complete  set  of  comparisons 
given  by  Rimer,  et  al.  (1979),  indicate  that  the  main  features  of  the 
ground  motion  have  been  modeled  rather  well.  The  main  discrepancy 
is  that  spallation  effects  are  too  large  within  300  meters  of  ground 
zero.  There  are  other  details  that  are  not  matched  by  the  calcula¬ 
tion,  but  much  of  the  discrepancy  can  reasonably  be  attributed  to 
the  elementary  axisymmetric  geometry.  Also,  little  information 
about  the  constitutive  properties  of  the  weathered  rock  near  the 
surface  was  available.  A  poor  model  for  this  layer  is  probably  the 
explanation  for  the  surface  spallation  being  too  large. 

Having  made  a  satisfactory  calculation  of  PILEDRIVER,  the 
yield  and  burial  depth  were  varied  and  three  more  calculations  were 
done.  These  are  the  four  S-Cubed  granite  calculations  to  be 
discussed  in  this  report. 

4.2  SOURCE  DESCRIPTION 

The  four  S-Cubed  granite  calculations  were  for  three  depths 
and  three  yields.  These  are: 


Depth 

Yield 

:4 

Identifier 

(m) 

(kt) 

P01 

460 

60 

• 

GRAN1 

1000 

150 

(  - 

GRAN2 

1000 

20 

4 

GRAN3 

400 

20 

The  source  region 

was  a  layered  half space 

meant  to 

represent 

• 

the  PILEDRIVER  environment.  This  structure  is 

given  in 

Table  6. 

Note  that  the  granite 

at  the  source  depths  is 

a  higher 

velocity 

material  than  the  ATI  granite  (a  *  4.4  km/sec. 

8  *  2.54 

km/sec) . 

The  shear  modulus  is  20cj 

larger  (206  kbar  compared 

to  172  kbar). 

> 
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tfi-SL 

X  VELOCITY 


Figure  23.  Measured  and  calculated  (heavy  line)  horizontal  velocities  at  shot  level 
station  16-SL  (range  =  470  ra).  The  data  marked  UH  are  from  a  velocity 
gauge  while  AH  indicates  an  integrated  accelerometer  record. 
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tegrating  an  accelerometer  record. 


The  velocities,  u'J,  ii^  (r,z,t),  and  stresses, 
cr'rr»  ^zz  ’  arz  (r»2»t)i  were  monitored  at  stations  on 
a  cylindrical  surface  surrounding  the  region  of  inelastic  material 
response.  The  radius  and  depth  to  the  bottom  of  this  cylindrical 
surface  are  listed  in  Table  7,  together  with  the  number  of 
monitoring  stations  on  the  side  and  bottom.  The  station  spacing  is 
about  25  meters  for  all  but  GRANl,  where  it  is  50  meters.  This 
compares  to  the  ATI  calculations  where  the  station  spacing  was  about 
50  meters. 

In  Figure  25  we  plot  some  representative  velocity  and  stress 
time  histories  for  each  of  the  calculations.  The  GRAN2  calculation 
is  very  close  to  completion,  since  the  velocities  are  nearly  zero  on 
the  monitoring  surface.  The  others  have  been  terminated  at  a  time 
when  motion  is  not  entirely  stopped  at  the  monitoring  surface.  The 
GRANl  is  probably  the  worst  in  that  resoect. 

Note  that  the  monitored  stresses  attain  static  values  that  are 
a  significant  fraction  of  the  peak  values.  This  is  in  contrast  to 
the  ATI  calculations  where  the  static  tractions  were  very  small. 
The  reason  is  the  much  smaller  scaled  distance  to  the  monitoring 
surface  in  the  S-Cubed  calculations  (the  radius  of  the  monitoring 
surface  in  the  ATI  calculations  was  about  2.1  kilometers).  The 
static  stresses  decay  approximately  as  1/R^. 

In  Figure  26  we  plot  the  values  of  the  monitored  stresses  at 
the  last  time  point  versus  position  on  the  monitoring  surface. 
These  are  remarkably  smooth  except  near  the  free  surface. 

4.3  THE  VERTICAL  e0RC5  AND  IMPULSE 

For  these  calculations  the  vertical  force  and  impulse  are 
given  by 


Fz(t) 


,/f, 


4<ro. 


b) 


r0dr0 


*  2. 


/ 


a  azr(a,z)  dz,  (20) 
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TABLE  6 


SOURCE 

REGION  STRUCTURE 

FOR  S-CUBED 

GRANITE 

CALCULATIONS 

Depth 

(km) 

Thickness 

(km) 

a 

(km/sec) 

8 

(km /sec) 

0 

(km/sec) 

0.05 

0.05 

1.44 

0.752 

2.65 

0.18 

0.13 

4.69 

2.46 

2.65 

00 

oo 

5.35 

2.79 

2.65 

TABLE  7 

MONITORING  STATIONS  FOR  S -CUBED  GRANITE  CALCULATIONS 


Source 

Monitoring  Surface 

Number  of 

Depth 

vield 

Depth 

Radius 

Stations 

Identifier  (km) 

(kt) 

(m) 

(m) 

Bottom  Side 

P01  460 

60 

1057 

1200 

45 

54 

>D 

O 

o 

o 

150 

1696 

1088 

22 

36 

GRAN2  1000 

20 

1391 

1088 

44 

61 

GRAN3  400 

20 

838 

1088 

44 

39 
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Figure  25.  The  monitored  velocities  and  tractions  are  plotted  at  selected 
positions  (r  and  z  in  kilometers).  For  consistent  sign  conven¬ 
tions,  the  u§  and  oPz  for  PD1  must  be  inverted. 
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Figure  25.  (continued) 
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Figure  2F.  The  monitored  stresses  at  the  last  time  point  are  plotted 
versus  position.  The  abscissa  i$  actually  the  node  num¬ 
ber,  starting  with  1  at  the  axis  of  symmetry  and  counting 
across  th<»  bottom  ana  up  the  siae  of  the  cylirarical 
surface. 
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where  a  is  the  radius  and  b  is  the  depth  of  the  cylindrical 
monitoring  surface.  The  analogous  equations  for  the  ATI  calcu¬ 
lations  are  (1)  in  Section  2.3.  In  Sections  2.3  and  2. a  we  discuss 
the  downward  force  and  impulse  in  a  general  way  and  point  out  that 
conservation  of  linear  momentum  requires  that  they  vanish  at  late 
time.  A  procedure  for  applying  a  correction  factor  to  the  vertical 
tractions  to  force  the  calculation  to  satisfy  these  conditions  was 
outlined  in  Section  2. a.  in  this  case,  we  have  the  stresses  rather 
than  the  tractions,  and  the  implementation  is  slightly  different. 
Once  again,  we  assume  that 

FZ(T1)  -  K  , 

IZ(T1)  =  -L  ,  (21) 

where  t  *  Tj  at  the  last  time  step.  Then  if  K  and  L  have  the  same 
sign,  we  can  zero  the  force  and  impulse  in  the  following  way.  We 
add  another  time  step  at  =  T^  -  Tp  such  that 

Fz(T2)  =  Iz(T2)  =  0  *  {22) 

This  is  done  by  setting 

4!r-b’T2>  *  <’22,r'5'Tl>  *  '  °£ri3- 

v*(«.z.y  -  .rz<»,i.Tj) » ,£  ,  o<2<b,  <23) 

and  taking 


At  *  2L/K. 


(?*) 
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The  correction  terms  are  computed  from 


.  M  -K  I  M 

6a  *  ■  a 

zz  r  |  zz 


(25a) 


for  stations  on  the  bottom  of  the  cylindrical  monitoring  surface,  and 


,  M  -K 
5arz  =  T" 


rz 


(25b) 


for  stations  on  the  side,  where 


Wb'Ti>i  Vo  * 


a21“(a  ’  2*Tl ) 


adz  .  (26) 


The  other  quantities  are  extended  to  7^  as  follows: 


V  *  • 

ur(T2)  *  uz(‘l)  ’ 

£<V  ■  v<V  • 

(27) 

-  ’"z'v  • 

0  £  r  <  a  , 

■Vz'V  ■  • 

0  <  z  <  b  . 

If  K  and  L  have  opposite  signs,  it  is  mucn  more  difficult  to 
zero  the  impulse.  This  happens  to  be  the  case  for  three  of  the  four 
S-Cubed  calculations.  For  the  other,  the  at  from  (24)  is  too  large 
for  this  method  to  be  applied  very  satisfactorily. 

In  Figure  27  we  show  the  force  and  impulse  for  each  of  the 
S-Cubed  calculations.  The  force  is  nearly  zero  in  every  case,  but 
the  impulse  is  not.  This  is  because  there  is  a  substantial  amount 
of  momentum  remaining  inside  the  monitoring  surface  at  the  end  of 
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Figure  27.  The  vertical  force  and  impulse  for  the  four  S-Cubed  granite 
calculations  are  dotted.  A  dashed  line  indicates  the  imposed 
final  time  steo  wnich  zeroes  the  force  (the  line  is  more 
properly  plotted  as  a  Quadratic  for  the  imoulse). 
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the  calculation.  The  actual  momentum  within  the  cylindrical  volume 
was  calculated  and,  indeed,  is  in  good  (within  a  factor  of  1.5) 
agreement  with  the  IZ(T^)  in  Table  3,  which  was  calculated  from 
the  surface  integral  (20).  Most  of  the  momentum  is  in  material  that 
has  been  spalled.  This  can  be  seen  in  crack  density  plots  which 
will  be  discussed  later,  in  Section  5.7.  For  GRANl,  which  we  know 
(Figure  25)  is  furthest  from  completion,  the  dominant  momentum  is 
still  associated  with  upward  motion.  For  the  others,  the  positive 
impulse  at  the  final  time  indicates  that  the  spalled  material  is 
mainly  moving  downward,  closing  cracks. 

As  well  as  the  impulse  being  far  from  zero,  the  condition  that 
K  *  L  >  0  is  not  satisfied  in  two  of  the  four  cases.  Therefore,  we 
decided  to  allow 

Iz(T2)  *  0  .  (28) 

Rather  than  using  (24)  to  find  the  time  step,  we  pre-selected  at  to 
be  0.5  seconds  when  K  *  L  >  0  and  0.05  seconds  when  K  *  L  <  0. 

The  parameters  for  correcting  the  force  and  impulse  are  given 
in  Table  8.  The  corrections  to  the  stresses,  K/r  in  eauation  (25), 
are  not  very  large,  especially  for  PD1  and  GRAN2. 

Tne  inability  to  zero  the  impulse  means  that  the  monitored 
solutions  do  not  satisfy  conservation  of  linear  momentum.  As  is 
shown  in  Appendix  A,  this  means  that  the  contribution  of  the 
vertical  force  terms  is  larger  than  it  should  be  at  long  periods. 
However,  as  we  will  demonstrate  in  the  next  section,  the  error  turns 
out  to  be  insignificant  at  the  periods  of  interest. 
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TABLE  8 


PARAMETERS  FOR  CORRECTING  THE  VERTICAL  FORCE  AND  IMPULSE  TO  ZERO* 


w  w  w 

1 1  /  i  n!5  /  i  r\l  5 


Iden¬ 

tifier 

Depth 

(m) 

Yield 

(kt) 

T1 

(sec) 

(1015 

dyne- 

sec) 

(1015 

dyne- 

sec) 

(1015 

dyne- 

sec) 

(1016 

dynes) 

K/r 

at 

(sec) 

PD1 

460 

60 

1.68 

2.9 

14.6 

14.7 

10.0 

-0.29 

0.05 

GRANl 

1000 

150 

1.92 

33.6 

-23.5 

-15.1 

10.3 

0.31 

0.5 

GRAN2 

1000 

20 

1.68 

-5.0 

3.7 

2.4 

3.3 

-0.15 

0.5 

GRAN3 

400 

20 

1.43 

9.1 

9.7 

9.9 

3.2 

0.29 

0.05 

The  true  values  of  F2,  I2  and  are  obtained  by  multiplying  the  listed 
values  by  2tt. 
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V.  M  AND  m,  ESTIMATES  FOR  S-CUBED  GRANITE  CALCULATIONS 
S  D 

5.1  INTRODUCTION 

The  granite  source  calculations  described  in  the  orevious 
section  were  processed  with  the  methods  described  in  Appendix  A  to 
compute  synthetic  body  and  surface  wave  seismograms  in  the  same 
earth  models  used  with  the  ATI  sources  in  Sectton  III.  In  this  sec¬ 
tion  we  describe  these  calculations  and  analyze  the  results.  Tne 
ATI  and  S-Cubed  calculations  will  be  compared  in  Section  VII. 

We  begin  by  discussing  the  far-field  Rayleigh  waves  in  a  half¬ 
space.  These  are  easily  compared  to  the  simole  Rayleigh  pulse  we 
expect  for  an  explosion  source.  Analysis  of  these  Rayleigh  waves 
leads  to  the  conclusion  that  the  contribution  of  the  vertical 
tractions,  which  we  know  to  include  some  error,  is  not  very  large. 

Also,  we  use  the  halfsoace  Rayleigh  waves  to  determine  an  equivalent 
RDP  representation  for  the  two-dimensional  source.  Finally,  we 
compute  Rayleigh  waves  in  the  central  United  States  crustal  model 
used  earlier  with  the  ATI  sources. 

Tne  far-field  body  waves  are  shorter  period  and  are  associated 
with  wave  propagation  along  a  particular  takeoff  angle.  Thus,  they 
are  somewhat  easier  to  understand.  The  effects  of  the  nonlinear 
interaction  with  the  free  surface  can  be  directly  seen.  In  connec¬ 
tion  with  these  seismograms,  we  again  compute  an  equivalent  RDP  to 
directly  display  the  two-dimensional  effects.  We  also  look  at 
spallation  in  some  detail,  presenting  Dlots  of  the  crack  distri¬ 
bution  for  each  calculation. 

An  especially  important  feature  of  these  granite  calculations 
is  that  we  are  able  to  compare  with  one-dimensional  source  calcu¬ 
lations  done  with  the  same  constitutive  model.  It  is  satisfyina  to 
see  that  the  two-dimensional  effects  are  in  accordance  with  our 
intuitive  expectations,  being  rather  small  for  the  deep  events  and 
important  for  the  shallow  ones.  The  direct  (downward  propagating)  P 
wave  from  all  the  two-dimensional  calculations  is  about  the  same  as 
predicted  with  the  one-dimensional  source.  These  conclusions  are 
summarized  in  Section  5.10. 
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5.2  FAR— F I ELl)  RAYLEIGH  WAVES  IN  A  HALFSPACE 


As  with  the  ATI  calculations  (Section  3.2),  we  first  look  at 
the  Rayleigh  pulse  in  a  halfspace.  In  this  case,  the  source  calcu¬ 
lations  were  done  in  a  three-layered  model.  For  the  propagation 
medium  we  chose  a  halfspace  with  a  a  6.0  km/sec,  8  =  3.5  km/sec,  o  » 
2.7  gm/cm3.  The  techniaue  of  Bache,  Rodi  and  Harkrider  (1978), 
which  is  described  in  Section  A. 4,  was  used  for  the  transition 
between  the  two  structures.  The  range  was  1000  km  and  the  WWSSN  15 
-  100  instrument  response  was  included.  For  attenuation,  we  had  0  * 
10  in  the  top  kilometer  and  Q  =  300  elsewhere  in  the  propagation 
medium. 

In  Figure  28  we  compare  the  halfspace  Rayleigh  pulse  from  the 
S-Cubed  two-dimensional  calculations  to  the  Rayleigh  pulses  from  ROP 
sources  at  the  same  depths  and  yields.  The  RDP  sources  were 
computed  at  two  depths  with  the  same  constitutive  model  used  in  the 
two-dimensional  calculations.  The  overburden  pressure  difference 
between  400  and  460  meters  was  ignored.  These  RDP  source  functions 
are  plotted  in  Figure  29. 

The  Rayleigh  pulses  for  the  two-dimensional  calculations  have 
waveforms  much  like  those  from  the  ROP  sources.  There  is  no  obvious 
indication  of  numerical  difficulty.  In  the  next  two  sections,  we 
will  examine  this  comparison  more  closely. 

5.3  CONTRIBUTION  OF  VERTICAL  TRACTIONS 

As  pointed  out  at  the  end  of  Section  IV,  these  calculations 
could  not  easily  be  corrected  to  bring  the  total  vertical  force  and 
impulse  to  zero  at  late  time.  The  asymptotic  behavior  of  the 
solution,  discussed  in  Appendix  A,  indicates  that  the  contribution 
of  the  vertical  traction  terms  (arz  and  azz)  is,  therefore,  too 
large  at  long  period.  Thus  we  know  our  surface  wave  solution  is 
incorrect  at  very  long  periods,  but  the  importance  of  this  static 
offset  error  at  the  periods  of  interest  for  M  can  only  be  esti¬ 
mated  by  examining  the  numerical  results. 
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Figure  29.  The  spectral  amplitude  of  the  reduced  velocity  potential 
is  plotted  for  two  spherically  symmetric  source  calcula¬ 
tions  in  the  PILEDRIVER  source  material.  The  calcula¬ 
tions  differ  only  in  the  overburden  pressure.  The  ampli¬ 
tude  axis  is  scaled  to  0.02  KT,  while  the  frequency  axis 
is  scaled  to  60  KT. 
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igure  28.  ^alfsoace  Rayleign  waves  for  the  S-Cubed  calculations  are  com¬ 
pared  to  tnose  computed  with  an  POP  source.  Listed  with  eacn 
record  is  the  oeak-to-oeak  amplitude,  wnich  has  been  corrected 
for  the  instrument  response  at  the  aooarent  oeriod  IT). 
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In  Figure  30  we  comoare  the  halfspace  Rayleigh  wave  solution 
from  Figure  28  to  the  Rayleigh  wave  solution  ccmDuted  from  only  the 
vertical  traction  contribution.  The  Rayleigh  wave  spectra  for  these 
seismograms  are  plotted  in  Figure  31.  The  period  of  the  Rayleigh 
pulses  in  Figure  30  is  5  to  8  seconds,  and  the  spectral  amplitude 
differences  at  those  periods  are  fairly  represented  by  the  peak- 
to-peak  amplitudes  of  the  seismograms. 

On  each  of  the  spectral  plots  in  Figure  31,  we  have  marked  the 
frequency  associated  with  the  time  (T^)  at  which  the  calculation 
was  terminated.  As  a  rough  approximation,  we  can  say  that  the 
solution  for  frequencies  higher  than  this  is  as  accurate  as  the 
finite  difference  calculations,  while  moving  to  lower  frequencies 
tends  to  exaggerrate  the  effect  of  errors  in  the  calculated  values 
at  the  last  time  step. 

The  contribution  of  the  vertical  force  term  is  small  at  the 
periods  of  interest  for  the  two  deep  source  calculations,  GRANl  and 
GRAN2.  This  is  consistent  with  the  ATI  calculations,  Figure  8, 
where  the  vertical  force  contribution  was  relatively  more  important 
at  shallow  depths,  even  when  the  impulse  had  no  static  offset.  Only 
for  GRAN3  does  the  (incorrect)  static  offset  in  the  vertical  impulse 
appear  to  be  causing  significant  errors.  But  even  in  this  case,  the 
Rayleigh  wave  solution  is  probably  accurate  for  periods  out  to  20 
seconds  or  so. 

5. A  AN  ECU  I VALENT  ROP  FROM  SURFACE  WAVES 

As  discussed  in  appendix  3  and  in  Section  3. a,  another 
interesting  way  to  display  the  propagated  surface  waves  is  in  terms 
of  an  equivalent  ROP.  That  is,  we  comoute  the  ROP  source  which 
gives  the  same  Rayleigh  waves  as  the  two-dimensional  source,  ^rom 
(16)  in  Section  3. A,  we  see  that  the  equivalent  RDP,  fg,  is  simply 
the  ratio  of  the  spectra  of  the  seismograms  in  Figure  28  times  the 
RDP  used  for  those  calculations,  which  are  plotted  in  Figure  29. 
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Figure  30.  The  complete  halfspace  Rayleigh  waves  for  the  S-Cubed  granite 
calculations  are  compared  to  the  solution  comouted  with  only 
the  vertical  traction  terms.  The  peak-to-oeak  amplitudes  in 
microns  at  15  seconds  are  shown  with  each  record. 
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PILEDRI VER 
W  =  60  KT 
M  =  0.46  kin 


Figure  31.  The  Rayleigh  wave  spectra  are  plotted  for  the  seismograms  of  Figure  30.  These  spectra  do  not 
include  the  effect  of  Q  or  the  seismometer.  The  total  spectrum  is  the  larger  at  long  period 
on  each  plot. 


In  Figure  32  we  plot  the  if  l  for  the  four  S-Cubed  granite 
calculations.  They  are  compared  to  l*!  for  the  RPP  sources  *69  or 
472,  scaled  to  the  appropriate  yield.  Our  remarks  about  the  high 
freouency  behavior  of  l1*1  |  in  Section  3.4  apply  here,  again,  the 
spectral  "holes"  shift  to  lower  frequency  with  depth  and  are 
believed  to  be  associated  with  the  nonlinear  interaction  with  the 
free  surface,  including  spallation.  We  also  point  out  that  the  long 

1 4  | 

period  errors  expected  in  GRAN3  are  apparent  in  the  i  (it 
increases  for  periods  longer  than  about  20  seconds).  For  the  other 
calculations,  the  to  !  are  nearly  constant  at  long  period,  as 
they  should  be.  This  is  another  reason  to  believe  that  our  failure 
to  zero  the  vertical  impulse  is  not  important  at  the  oeriods.  of 
interest. 

!  A  i 

Another  display  of  thei'i'g1  for  the  S-Cubed  granite  calcula¬ 
tions  is  shown  in  Figure  33.  The  1$  |  for  all  four  calculations 
are  plotted  together  at  the  top  of  the  figure.  It  is  the  long 
period  behavior  (called^)  that  is  of  most  interest.  For  the  two 
deep  explosions  the  ^  deduced  from  the  2-0  calculations  is  about 
the  same.  Further,  this  4^  is  nearly  the  same  as  that  obtained  from 
a  1-D  calculation  at  the  same  depth.  This  is  an  important  result 
that  gives  considerable  confidence  in  the  consistency  of  the  whole 
procedure.  If  we  refer  back  to  the  half space  Rayleigh  wave 
comparison  in  Figure  28,  we  also  see,  as  expected,  close  agreement 
in  the  amplitudes  of  the  RDP  and  2-D  seismograms  for  the  deep 
events,  GRAN1  and  GRAN2.  At  shallow  depths,  the  two-dimensional 
calculation  predicts  much  larger  (a  factor  of  2  to  3)  surface  waves 
than  a  one-dimensional  calculation.  We  exoect  2-D  effects  to  be 
more  important  for  shallow  source  depths,  so  this  is  not  surorising. 

At  the  bottom  of  Figure  33  the  1-0  and  2-0  source  functions 
for  P ILEDR TVER  are  compared  to  the  Mueller /Murphy  granite  source  for 
PILEDRIVER.  The  S-Cubed  1-0  source  has  nearly  the  same  'V  ,  but  is 

-30 

very  different  at  frequencies  above  0.1  Hz. 
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Figure  32.  The  j'i'gj  for  the  S-Cubed  granite  calculations  are  comoared  to 
source  469  or  472  (Figure  29),  deoending  on  the  deoth.  The 
amplitudes  have  been  scaled  to  a  common  yield,  0.02  KT  and  the 
frequencies  have  been  scaled  to  60  KT. 
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Figure  33.  The 


The  spectral  amDlitude  (scaled  to  0.02  Kt)  of  the  equivalent 
R VP,  |ts|,  is  plotted  for  each  of  the  S-Cubed  aranite  calcu¬ 
lations.  At  the  bottom  the  p ILEDRI VER  |?e  |  is  comoared  to 
1-0  source  469  and  the  Muel ler /Murphy  granite  source  at  this 
yield.  The  frequency  axes  are  scaled  to  60  kt. 
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5.5  FAR— FIELD  RAYLEIGH  WAVES  IN  REALISTIC  EARTH  MODELS 

The  fundamental  mode  Rayleigh  waves  were  comouted  for  the  four 
S-Cubed  granite  calculations  with  the  same  path  model  used  for  the 
ATI  sources  (Section  3.5).  This  model  is  listed  in  Table  2,  with 
the  exception  that  the  top  three  kilometers  of  the  source  region 
model  were  replaced  by  the  source  structure  in  Table  6. 

The  computed  Rayleigh  waves  are  shown  in  Figure  34.  These  are 
directly  comparable  to  the  surface  wave  synthetics  in  Figure  11  for 
the  ATI  source  calculations.  The  values  were  computed  ac¬ 
cording  to  (13)  and  are  listed  in  Table  9.  These  values  will  be 
discussed  in  more  detail  in  Section  VII  where  we  compare  them  with 
the  from  the  ATI  sources. 

5.6  FAR-FIELD  30DY  WAVES 

Short  period  body  wave  seismograms  like  those  for  the  ATI 
sources  in  Figure  15  are  plotted  in  Figure  35.  We  also  plot  the 
analogous  seismograms  for  the  1-D  sources  469  (PILEDRIVER  and  GRAN3) 
and  472.  The  only  difference  between  the  calculations  here  and 
those  for  the  ATI  sources  in  Section  3.6  is  that  the  top  three 
kilometers  of  the  source  region  structure  (Table  4)  are  replaced 
with  the  structure  in  Table  6. 

The  body  wave  magnitude  values  for  the  calculations  in  Figure 
35  are  listed  in  Table  10.  These  m^  were  computed  using  (19). 
Note  that  the  magnitudes  for  the  1-D  source  are  very  nearly  the  same 
as  those  for  the  2-D  sources,  especially  for  the  two  deep  events. 

For  the  two  deeo  sources,  GRANl  and  GRAN2,  the  seismograms  for 
the  ROP  and  two-dimensional  sources  have  very  similar  waveforms. 
The  two-dimensional  effects  are  much  more  prominent  on  the  seismo¬ 
grams  for  the  shallow  sources.  This  is,  of  course,  what  we  should 
expect. 
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he  S-Cubed  granite  calculations. 
seismoaieter  response  is  included. 


TABLE  9 


Ms  FOR  S-CUBED  GRANITE  CALCULATIONS 


Calculation 

Depth 

(km) 

Yield 

C<t) 

A 

(microns ) 

T 

(sec) 

Ms 

PILEDRIVER 

0.46 

60 

8.7 

13  3 

4.55 

GRANl 

1.0 

150 

8.7 

12.4 

4.44 

GRAN2 

1.0 

20 

1.0 

13.0 

3.56 

GRAN3 

0.40 

20 

3.3 

13.9 

4.14 
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Figure  35.  Body  wave  seismograms  for  the  S-Cubed  granite  calculations 
are  compared  to  seismograms  with  an  RDP  source.  The  peak- 
to-peak  amplitudes  listed  with  each  seismogram  have  been 
corrected  for  the  instrument  response  at  the  apparent 
period  (T  )  in  Table  10. 
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TABLE  10 


mb  FOR  S-CUBED  GRANITE  CALCULATIONS 


Calculation 

Depth 

(km) 

Yield 

(kt) 

b 

(nm) 

T. 

0 

(sec) 

b 

mb 

C 

(nm) 

TC 

(sec) 

■s 

PILEDRIVER 

0.A6 

60 

TWO-O IMF MS IONAL 

169  0.73  5.62 

130 

0.75 

5.63 

GRANl 

1.0 

150 

553 

0.81 

6.09 

936 

0.85 

6.29 

GRAN2 

1.0 

20 

57 

0.67 

5.18 

33 

0.70 

5.32 

GRAN3 

0.4 

20 

64 

0.71 

6.21 

120 

1.07 

5.30 

ONE-DIMENSIONAL 


469 

0.46 

60 

166 

0.67 

5.65 

247 

0.70 

5.30 

472 

1.0 

150 

568 

0.73 

6.11 

999 

0.82 

6.33 

472 

1.0 

20 

52 

0.68 

5.22 

98 

0.72 

5.3° 

469 

0.4 

20 

44 

0.60 

5.11 

58 

0.65 

6.21 
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5.7  ANALYSIS  OF  SPALLATION 


Probably  the  most  important  two-dimensional  effect  is  the 
occurrence  of  spallation.  This  complex  phenomenon  degrades  and 
distorts  the  free  surface  reflected  waves  and  partitions  some  of  the 
energy  into  spall  closure  phases.  This  is  an  extremely  complex 
process  which  can  be  seen  by  examining  displays  of  the  crack 
porosity.  We  will  show  some  attempts  to  reproduce  the  synthetic  for 
the  2-D  source  by  modifying  the  seismogram  for  the  1-D  source  by  (1) 
suppressing  pP  by  a  constant  factor,  and  (2)  adding  an  impulsive 
force  at  the  surface  to  represent  spall  closure.  This  was  not  very 
successful,  apparently  because  this  model  is  much  too  simole  to 
represent  this  calculation. 

In  Figure  36  we  show  crack  density  plots  for  the  PILEDRIVER 
calculation  at  three  times.  These  and  similar  plots  for  the  other 
three  calculations  (Figures  37  through  39)  are  in  two  different 
formats  which  are  explained  in  the  caption.  The  plot  at  0.426 
seconds  shows  the  maximum  extent  of  cracking  near  the  surface.  The 
entire  first  layer  (Table  5),  which  is  50  meters  thick,  is  "spalled" 
to  a  radius  of  1  kilometer.  The  spalled  region  would  have  extended 
to  even  greater  radii,  but  was  terminated  at  this  radius  to  ensure 
elastic  behavior  at  the  monitoring  surface.  Recall  that  in  Section 
4,1  we  mentioned  that  the  amount  of  spall  is  unreal istically  large 
in  the  calculation,  probably  due  to  an  inaccurate  material  model  for 
the  surface  layer. 

The  plots  at  0.63  seconds  show  that  some  of  the  near  surface 
cracks  have  closed.  Spall  closure  pulses  corresponding  to  these 
cracks  can  be  seen  on  the  records  at  the  free  surface  monitoring 
stations  (e.g..  Figure  24).  At  1.63  seconds,  the  end  of  the  cal¬ 
culation,  many  more  of  the  near  surface  cracks  have  closed,  but 
there  is  a  substantial  region  still  in  free  fall.  Cracked  material 
is  distributed  extensively  through  the  grid  even  at  this  time.  It 
is  clear  that  spall  closure  is  very  complex  and  is  distributed  over 
more  than  a  second  in  time.  Certainly  the  pP  must  be  degraded  to  a 
substantial  degree  in  opening  these  cracks. 
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Figure  36.  The  crack  density  is  plotted  at  three  times  for  the  P ILEDRI VER 
calculation.  Dots  or  open  boxes  indicated  hoop  (out-of-plane) 
cracks.  In-plane  cracks  open  across  the  lines.  The  X  or 
closed  boxes  indicate  zones  in  which  tensile  failure  has 
occurred  in  all  three  coordinate  directions.  These  zones  are 
thus  decoupled  from  adjacent  zones. 
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Cycle:  900 


Time  (msec)  566 


38.  The  crack  density  is  plotted  at  two  times  for 
GRANT  (VI  *  150  KT  and  H  =>  T.O  km). 
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Cycle:  2500 


Time  (msec)  1887 


Figure  38.  (continued) 
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In  Figure  40,  we  again  compare  the  seismograms  for  the  one- 
and  two-dimensional  PILEDRIVER  sources  (Figure  35),  then  repeat  the 
comparison  with  the  upgoing  waves  entirely  suppressed.  We  see  that 
suppressing  the  pP  does  improve  the  agreement,  especially  of  the 
relative  amplitudes  of  the  three  main  peaks.  The  bottom  two 
comparisons  in  the  figure  will  be  discussed  in  later  paragraphs. 

In  Figure  37  we  show  the  crack  density  at  several  times  for 
GRAN2,  which  was  an  overburied  explosion.  This  event  had  the  least 
cracking  of  the  four.  Even  so,  we  see  that  the  50  meter  thick  top 
layer  is  entirely  spalled.  Some  cracks  remain  open  at  the  last 
time,  but  not  nearly  as  many  as  for  P ILEDR I VER . 

In-  Figures  38  and  39  we  show  crack  density  plots  for  GRANl  and 
GRAN3.  Again,  the  first  layer  is  entirely  spalled.  Many  cracks 
remain  open  at  the  end  of  the  calculation,  and  the  cracked  region 
extends  to  considerable  depth.  These  two  calculations  are  between 
PIIEDRIVER  and  GRAN2  in  tne  extent  of  cracking. 

We  should  remark  that  the  last  crack  density  plot  for  each 

calculation  shows  the  cracks  that  remain  open  when  the  calculation 
is  terminated.  There  is  clearly  a  substantial  amount  of  momentum 
remaining  inside  the  monitoring  surface  at  this  time.  Thjis 
probably  accounts  for  most  of  the  "static"  offset  of  the  vertical 
impulse  discussed  in  Section  4.3. 

How  much  seismic  energy  is  generated  by  the  cracks  that  do 

close?  In  all  four  calculations,  the  entire  first  layer  is  spalled 
to  a  radius  of  1000  meters.  In  most  cases  the  spalled  region  is 

substantially  thicker  in  the  region  directly  above  the  explosion. 
Tne  spall  closure  phases  from  such  a  large  mass  of  material  would 
make  a  huge  seismic  signal  if  they  were  in  phase. 

The  volume  of  a  cylinder  with  radius  1000  m  and  thickness  50  m 

12  3  3 

is  50*  x  10  cm  .  The  density  is  2.65  gm/cm  ,  so  the  mass  of 

this  volume  is 

M  »  4.16  x  1014  gm.  (29) 
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Figure  40.  The  seismogram  for  the  two-dimensional  PIIEDRIVER  cal¬ 
culation  is  compared  to  synthetic  seismograms  constructed 
v*i th  elementary  point  sources.  The  parameters  for  Ps  are 
1.4  x  10^6  dyne-sec  for  the  amplitude  and  1.7  seconds  for 
the  lag  time. 
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If  we  follow  Viecelli  (1973)  and  assume  the  soall  seoaration  is 
linearly  proportional  to  range,  the  total  impulse  is 

I=§(2gh0)1/2  ,  (30) 

where  g  is  the  gravitational  acceleration  and  hQ  is  the  maximum 
height  of  spall.  But  if  hQ  is  only  20  cm,  we  get 

I  =  2.75  x  1016  dyne-sec,  (31) 

which  happens  to  be  Viecelli's  estimate  for  the  spall  impulse  for  an 
event  of  this  yield  at  Pahute  Mesa.  The  actual  displacement  of  the 
spalled  zones  is  several  meters  (Rimer,  et  aj_,  1979). 

In  the  bottom  two  seismogram  pairs  in  Figure  40  we  compare  the 
2-0  PILFDRIVER  calculation  to  synthetic  seismograms  computed  with 
the  RDP  source  469  in  which  pP  is  partially  suDpressed  and  an 
impulsive  point  source  has  been  added  at  the  surface  to  represent 
soall  closure.  The  amplitude  of  the  spall  impulse  is  about  half  the 
value  in  (31),  so  we  see  that  an  impulse  that  large  would  make  a 
substantial  contribution  to  the  seismogram.  The  waveform  comparison 
in  Figure  40  is  not  very  good,  indicating  that  the  impulsive  ooint 
source  model  for  soall  closure  is  too  simple. 

The  point  of  this  analysis  is  that  examination  of  the  dimen¬ 
sions  of  the  SDalled  region  and  the  amount  it  is  displaced  suaaests 
that  coherent  spall  closure  would  give  very  large  seismic  phases. 
However,  a  phase  that  appears  to  be  associated  with  spall  closure  is 
only  seen  on  the  shallow  source  (PILEDRIVER  and  GRAN3)  records  and 
it  is  not  very  large  there.  Thus,  the  potential  energy  in 
spalled  material  in  these  calculations  is  not  very  efficiently 
converted  to  seismic  waves. 


5.8  800Y  WAVE  SPECTRA 

In  Figure  35  we  comoared  body  wave  seismograms  for  one-  and 
two-dimensional  calculations  of  the  four  sources.  One  steo  in 
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synthesizing  these  seismograms  is  to  compute  the  far-field  disolace- 
ment  spectrum  at  the  appropriate  takeoff  angle  (20.3*)  at  the  base 
of  the  source  structure  (Table  5).  These  far-field  diSDlacement 
spectra  are  compared  in  Figure  41.  The  SDectra  for  the  ROP  source 
have  the  familiar  shape  that  occurs  when  pP  is  a  soectral  replica  of 
P.  The  spectral  peaks  and  holes  occur  at  multiDles  of  the 
reciprocal  of  the  P-pP  delay  time. 

As  we  examine  the  spectra  for  the  two-dimensional  sources,  we 
see  that  the  spectral  scalloping  is  much  weaker  or  absent.  The 
overburied  event  GRAN2  shows  some  of  the  same  character  as  that  from 
the  analogous  RDP  source,  but  even  here  the  peaks  and  troughs  are 
less  prominent.  The  spectra  for  the  other  events  have  much  weaker 
and  less  regular  phase  interference  features.  These  spectra  are 
evidence  of  the  previously  stated  conclusion  that  the  pf  for  the  2-0 
calculations  is  relatively  small  and  is  not  a  spectral  replica  for 
°.  Spectral  interference  due  to  much  later  arriving  phases,  which 
we  associate  with  spall  closure,  is  also  seen  in  the  shallow  source 
spectra. 

5.9  AM  EQUIVALENT  RDP  FROM  BODY  WAVES 

As  with  the  surface  waves,  we  can  display  the  two-dimensional 
effects  in  a  direct  way  by  computing  the  RDP  that  gives  the  same 
seismogram.  Again,  we  use  (16)  in  Section  3.4.  That  is,  we  divide 
the  2-0  source  spectra  in  Figure  41  by  the  1-0  source  soectra  and 
multiply  by  the  RDP  source  function  used  for  the  latter. 

■Hie  spectral  amolitudes  of  the  eouivalent  R'/°  for  the  four 
S-Cubed  calculations  are  plotted  in  Figure  4?,  together  with  the 
aoorooriate  one-dimensional  source  functions.  The  high  frequency 
soikes  in  the  If!  are  due  to  the  presence  of  spectral  nulls  in  the 

1- D  source  soectra  that  are  not  replicated  in  the  2-0  spectra. 

Our  interest  in  body  waves  is  mainly  in  the  spectral 
amolitudes  of  the  source  function  near  1  Hz.  Note  that  the  1-0  and 

2- 0  source  functions  are  not  very  different  at  that  frequency.  We 
saw  this  in  a  different  way  when  comparing  time  domain  amolitudes  in 
Figure  35.  Also,  ignoring  the  soikes,  the  high  frequency  rolloff  of 
the  ;'f  is  about  that  of  the  analogous  1-0  source. 
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Figure  41.  P  wave  spectra  are  compared  for  the  one-dimensional  and  two- 
dimensional  source  calculations.  The  plots  are  in  log-log 
format  with  the  frequency  scaled  to  60  KT. 
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The  IV  i  estimated  from  body  waves  in  Figure  £2  should  be 
compared  to  estimates  of  the  same  quantity  from  surface  waves  in 
Figure  32.  Tnere  are  differences,  as  should  be  exDected.  if 
two-dimensional  effects  are  important,  the  jw  |  should  vary  with 
takeoff  angle  from  the  source. 

The  least  variation  between  body  and  surface  wave  estimates 
for  Jfg[  seems  to  be  for  the  two  deep  explosions,  GRAN1  and  GRAN?. 
Both  estimates  agree  that  two-dimensional  effects  are  not  very 
important  for  these  deep  events.  For  the  shallow  sources,  the 
equivalent  RDP  estimated  from  surface  waves  is  rather  different  from 
that  estimated  from  body  waves,  though  they  do  agree  that  the  2-0 
source  has  much  more  long  period  energy. 

5.10  SUMMARY 

Body  and  surface  wave  seismograms  computed  for  the  four 
S-Cubed  granite  calculations  give  consistent  and  believable 
results.  Our  main  conclusions  are: 

«  While  the  calculations  do  not  ideally  satisfy 
conservation  of  momentum  (the  total  vertical 
impulse  has  a  static  offset),  there  is  no 
indication  that  this  has  more  than  a  minor  effect 
on  the  solutions. 

•  Two-dimensional  effects  are  not  very  important 
for  the  two  deep  explosions.  Soth  m^  and 
are  little  different  from  the  values  estimated 
from  an  SOP  source  comouted  with  the  same  con¬ 
stitutive  model  at  the  same  death.  This  is  true 
even  though  considerable  cracking  and  soallation 
occur. 

•  For  the  shallow  events,  the  two-dimensional  ef¬ 
fects  act  to  enhance  the  excitation  of  surface 
waves  by  a  factor  of  two  or  three. 


127 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


•  For  body  waves,  the  P  phase,  which  controls 

m?,  is  essentially  the  same  for  one-  and 
two-dimensional  source  calculations. 

•  For  the  deep  events  both  m^  and  m£  are 

nearly  the  same  for  the  one-  and  two-dimensional 
sources.  However,  examination  of  spectra  in¬ 
dicate  that  pP  is  not  a  spectral  shadow  of  P. 
Tnis  is  not  surprising  in  view  of  the  cracking 
that  occurs  to  significant  aepths. 

•  For  shallow  events,  the  two-dimensional  effects 

substantially  alter  the  portion  of  the  waveform 
after  the  initial  P  wave  arrival.  Tnis  influ¬ 
ences  the  commonly  used  magnitude,  mp, 

though  by  less  than  0.2  m.  units. 
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VI.  COMPARISON  OF  PREDICTED  AND  OBSERVED  SEISMOGRAMS 

FOR  PILEDRIVER 

6.1  INTRODUCTION 

One  of  the  four  S-Cubed  granite  calculations  described  in 
Section  IV  was  intended  to  model  the  PILEDRIVER  event.  Rimer,  et 
al .  (1979)  showed  that  this  two-dimensional  calculation  gives  good 
agreement  with  the  main  features  of  the  recorded  near-field  ground 
motions.  In  this  section  we  will  make  a  careful  comparison  of  the 
far-field  body  and  surface  waves  from  this  calculation  with 
observations  of  PILEDRIVER. 

6.2  COMPARISON  OF  PREDICTED  AND  08SERVED  SURFACE  WAVES 

In  Section  5.5,  Figure  34,  we  showed  a  teleseismic  surface 
wave  for  the  PILEDRIVER  calculation.  The  M  was  comouted  from  a 
large  Airy  phase  that  is  a  feature  of  the  Darticular  structure  that 
was  used  for  the  synthesis.  This  M$  is  4.55,  compared  to  the 
observed  value  for  PILEDRIVER  of  about  3.88  (Eisenhauer,  unpublished 
work).  This  synthetic  M  value  is  actually  much  larger  than  we 
would  get  from  an  "average"  path  model,  as  we  will  demonstrate. 

Bache,  Rodi  and  Harkrider  (1978)  and  Bache,  Rodi  and  Mason 
(1973)  studied  the  Rayleigh  waves  from  NTS  explosions,  including 
PILEDRIVER,  at  the  WWSSN  stations  ALO  and  TUC.  In  the  former  study 
crustal  structures  were  inferred  for  the  NTS-ALO  and  NTS-TUC  paths. 

The  latter  study  was  concerned  with  estimating  the  source 
amplitudes  of  NTS  explosions,  using  the  crustal  models  to  correct 
for  the  path. 

In  an  earlier  report,  Bache,  Gouoillaud  and  Mason  (1977) 
compared  the  M  measured  by  Eisenhauer  with  an  M$  computed  from 
the  Airy  phase  amplitudes  at  ALO  and  TUC.  For  eighteen  events  at 
Yucca  Flat  and  °ahute  Mesa,  the  residual  was  minimized  by  using  the 
formulas: 

M^  ,  log  A  +  2.72  , 
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for  ALQ,  and 


ja 

I 


Mj  *  log  A  +  2.17  (32) 

for  TUC.  For  P ILEDRI VER ,  these  formulas  give 
,  3.82  , 

M l  .  3.65  ,  (33) 

which  are  reasonably  close  to  the  "mean"  value  of  3.38. 

Using  the  structures  from  Bache,  Rodi  and  Harkrider  (1978),  we 
computed  synthetic  seismograms  at  ALQ  and  TUC  using  the  2-0 
P ILEDR I VER  source.  The  results  are  compared  to  the  data  in  Figure 
43.  As  was  shown  in  Section  V,  the  2-0  effects  do  not  have  much 
influence  on  the  waveform.  The  Ms  values  computed  from  (32)  are 
4.19  for  ALO  and  3.93  for  TUC.  These  values  are  0.36  and  0.62  M 
units  smaller  than  value  (4.55)  we  computed  for  the  central  U.S. 
structure  in  Section  5.5. 

There  are  some  other  interesting  aspects  of  the  seismogram 
comparison  in  Figure  43  that  should  be  discussed.  The  synthetic  and 
observed  waveforms  are  in  excellent  agreement  at  ALO.  At  TUC,  the 
agreement  is  not  very  good,  even  though  the  synthetic  and  observed 
records  have  almost  the  same  phase  and  group  velocity  dispersion. 
This  poor  waveform  agreement  at  TUC  seems  to  be  due  to 
charactari sties  of  the  P ILEDR I VER  source.  cigure  44,  reproduced 
from  Bache,  Rodi  and  Harkrider  (1973),  compares  synthetic  and 
observed  '•ecords  for  typical  events  at  Yucca  "lat  and  Pahute  Mesa, 
as  well  as  for  PILEDR  TVER.  Tine  P  ILEDR  TVER  recording  at  TUC  is 
unique  in  its  disagreement  with  the  synthetic. 

Another  interesting  facet  of  the  comparison  in  Figure  43  is 
that  the  amplitude  ratio  is  not  very  different  at  the  two  stations. 
3ache,  Rodi  and  Mason  (1978)  found  that  the  observed  amplitudes  at 
ALO  were  relatively  small  compared  to  the  synthetics  for  an  ROP 
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Figure  43.  Observations  of  PILEORIVER  from  ALQ  and  TUC  are  compared  to 
synthetic  seismograms.  The  synthetics  are  on  the  bottom  for 
each  pair.  The  two  sets  at  the  too  of  the  page  are  from 
Bache,  Rodi  and  Harkrider  (197S)  and  were  computea  with  an  RDP 
source.  The  other  synthetics  were  computed  from  the  two-di¬ 
mensional  PILEORIVER  calculation.  The  peak-to-peak  amplitude 
in  microns  is  listed  with  each  record.  The  time  with  respect 
to  origin  times  for  the  ALQ  and  TUC  records,  observeo  and 
synthetic,  is  called  T0. 
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Figure  44.  Theoretical  and  observed  seismograms  are  compared  at  ALQ  (left)  and  TUC  for  events  in 
three  test  areas  at  NTS.  A  bar  indicating  one  minute  is  shown.  In  each  pair  the  ob¬ 
served  and  theoretical  records  start  at  the  same  time  with  respect  to  the  explosion 
detonation  and  this  time  is  indicated  as  T0.  The  theoretical  seismograms  were  com¬ 
puted  with  the  same  RDP  source  (reproduced  from  Bache,  Rodi  and  Harkrider,  1978). 


I 

source,  just  as  they  are  here.  But  in  that  study,  an  RDP  source 
that  matches  the  data  at  TUC  was  found  to  give  synthetic  seismograms 
that  are  50-60  percent  too  large  at  ALO  and  this  result  seemed 
independent  of  test  area.  However,  for  the  2-0  source,  the 
amplitude  ratio  difference  is  only  2.3/1. 9  or  20  percent. 

An  entirely  satisfactory  explanation  of  PILEDRIVER  surface 
waves  must  deal  with  the  observed  radiation  pattern,  which  has  been 
interpreted  as  evidence  for  a  double-couple  component  being 
superimposed  on  the  explosion  source  (Toksoz  and  Kehrer,  1972). 
Both  ALO  and  TUC  are  on  oositive  lobes  for  the  Toksoz  and  Kehrer 
(1972)  double-couple  fit  to  the  data,  but  at  azimuths  where  the 
inferred  double-couple  Rayleigh  waves  are  larger  than  those  from  the 
explosion  by  factors  of  3  to  3.5.  This  solution  suggests  that  the 
explosion  itself  couples  rather  weakly  into  surface  waves.  Such  a 
result  is  certainly  inconsistent  with  current  models  for  the 
constitutive  behavior  of  granite  represented  by  the  calculations  in 
this  report. 

Rivers  and  von  Seggern  (1979)  calculated  the  seismic  moment 
tensor  that  fits  the  observed  surface  waves  for  PILEDRIVER.  They 
concluded  that  the  orientation  of  the  superimposed  double-couple  is 
consistent  with  vertical  dip-slip  faulting,  rather  than  strike-slip 
faulting.  This  double-couple  reduces  the  observed  surface  waves  by 
about  a  factor  of  two.  The  two  estimates  for  the  PILEDRIVER 
double-couple  result  in  surface  wave  amplitudes  that  differ  by  a 
factor  of  six. 

If  the  Rivers  and  von  Seggern  (1979)  double-couple  solution 
were  added  to  our  two-dimensional  source  calculation  (with  little 
ohase  delay),  the  observed  and  synthetic  seismogram  amplitudes  in 
Figure  43  would  be  in  very  close  agreement.  Some  additional  source 
complexity  could  also  explain  the  waveform  differences  at  TUC. 

6.3  COMPARISON  OF  PREDICTED  AND  OBSERVED  BODY  WAVES 

We  begin  with  a  careful  comparison  of  synthetic  and  observed 
seismograms  at  two  particular  stations  that  we  have  used  in  previous 
synthetic  seismogram  studies.  We  then  compare  to  a  suite  of  WWSSN 
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data  collected  by  Hadley  and  Hart  (1979).  The  stations  to  be 
studied  in  detail  are  8FAK  (31uff,  Alaska,  a  *  33. A*)  and  HNME 
(Houlton,  Maine,  a  *  36.5*).  The  HNME  data  are  taken  from  the  LRSM 
shot  report  (SDL  Report  165)  for  this  event.  The  at  HNME 
is  5.53.  The  SDL  report  also  gives  a  mean  m^  of  5.43  for  stations 
with  a  >  16°.  At  8FAK  the  is  5.54.  These  values  are  not 
much  different  from  a  carefully  determined  network  average  mb  of 
5.47  given  by  Alewine,  Young,  Springer  and  Kleoinger  (unpublished 
report) . 

To  compute  synthetic  seismograms  at  specific  stations  we  must 
include  the  appropriate  upper  mantle  and  seismometer  responses.  For 
HNME  we  use  the  mantle  model  HWNE  (Helmberger  and  Wiggins,  1971)  and 
account  for  anelastic  attenuation  with  t*  *  0.8.  This  particular 
earth  model  and  t*  were  found  to  be  compatible  with  the  HNME  re¬ 
cordings  of  Pahute  Mesa  explosions  studied  by  Bache,  et  al .  (1979). 
At  this  range,  the  upper  mantle  response  is  not  too  different  from  a 

C  1 

constant  geometric  spreading  factor  of  6.4  x  10  km  ,  as  was 
mentioned  in  Section  3.6. 

For  BFAK  synthetics  we  use  the  path  model  selected  by  Bache, 
et  al .  (1975)  for  a  study  of  the  relative  body  wave  amplitudes  of 
NTS  explosions  in  different  test  areas.  This  upper  mantle  model  is 
a  slightly  modified  version  of  the  average  earth  model  C2  (Anderson 
and  Hart,  1976).  The  smooth  gradient  below  the  600  km  discontinuity 
was  modified  in  a  way  that  enhances  the  theoretical  amplitude  at  the 
range  of  3FAK  relative  to  the  original  C2  model.  Th is  amplitude 

response  is  roughly  eouivalent  to  a  constant  geometric  spreading 

-4-1 

factor  of  10  km  ,  which  is  25  percent  or  so  larger  than  the 
response  at  this  range  from  the  HNME  model  (3ache,  et  al.,  1975)  or 
the  Jeffrays-3ul len  average  earth  model  (Langston  and  Helmberger, 
1975). 

For  Alaskan  observations  of  NTS  explosions,  a  t*  of  1.05  was 
chosen  by  Bache,  et.  aJL_  (1975)  as  the  value  that  seemed  to  give  the 
best  agreement  with  the  observations  for  several  stations  and  many 
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events.  We  also  must  account  for  differences  in  the  seismometer 
response.  In  Figure  45  we  plot  the  response  for  3FAK  together  with 
the  LRSM  response  for  P ILEORI VER  at  HNME  and  the  KS360C0  reSDonse 
used  for  the  body  wave  synthetics  in  Sections  III  and  V. 

8efore  comparing  the  P ILEDRI VER  synthetics  to  the  HNME  and 
BFAK  data,  it  is  useful  to  look  at  how  this  comparison  turns  out  for 
some  other  events.  The  HNME  synthetic  seismograms  presented  by 
Bache,  et_  al  ♦  (1979)  were  for  three  representative  large  yield 

Pahute  Mesa  events,  MAST,  FONTINA  and  CAMEMBERT.  These  synthetics 
were  computed  with  the  Mueller /Murphy  tuff /rhyol ite  source  function 
(Mueller  and  Murphy,  1971),  modified  in  an  ad  hoc  way  to  account  for 
nonlinear  free  surface  interaction.  The  modifications  were  to 

multiply  the  upgoing  waves  from  the  source  by  a  constant  y  <  1 
(reducing  pP)  and  to  add  an  impulsive  downward  force  at  the  origin 
to  represent  spall  closure  (Table  11).  In  Section  5.6  we  showed 
some  seismograms  computed  in  this  way. 

In  Figure  46  and  Table  12  we  compare  synthetic  and  observed 
seismograms  at  HNME  and  BFAK  for  the  three  Pahute  Mesa  events.  The 
path  models  are  as  described  in  previous  paragraphs  and  the  source 
was  computed  with  the  y  and  spall  closure  phase  Ps>  described  in 

Table  11.  The  frequency  content  and  amplitudes  of  the  observations 
are  matched  very  well  at  HNME,  which  is  the  station  used  to  infer 
the  parameters  in  Table  11.  At  BFAK  the  freauency  content  is 
reproduced  pretty  well,  except  for  FONTINA  where  the  synthetic  has 
too  little  high  frequency  energy.  The  amplitudes  of  the  synthetics 
at  3FAK  are  too  large  by  40  oercent  or  so.  Considering  the  uncer¬ 
tainties  in  the  source  and  path  models,  these  are  auite  consistent 
results. 

In  Figure  47  and  Table  13  we  compare  the  3PAK  and  HNME 
observations  of  PILEDRIVER  with  synthetics  computed  with  our 
two-dimensional  source.  The  observations  for  this  event  are  in 

analog  form,  so  our  measurements  do  not  have  the  precision  they  did 
for  the  Pahute  Mesa  events.  At  BFAK  the  waveform  and  frequency 


PERIOD  (sec) 


Figure  45.  Seismometer  amplitude  response  for  the  Benioff  (LRSM)  and 
KS36000  seismometers  used  at  HNME  and  the  response  of  the 
SFAK  seismometer. 
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TABLE  11.  PARAMETERS  FOR  THE  Ps  PHASE  IN  THE 
SYNTHETIC  SEISMOGRAM  CALCULATIONS 


MAST 

CAMEMBERT 

F0NTINA 

Spall  Impulse 
( Dyne-sec /kt) 

10  x  1014W 

14  x  1014W 

9  x  1014W 

P-Ps  lag  time 
(sec) 

1.25 

1.92 

1.86 

Y* 

0.60 

0.50 

0.50 

Estimated  of  spall  impulse  from  near-field  data 
Viecelli  (1973):  3.5  x  1014W 

Sobel  (12978):  21-25  x  1014W 

Estimated  of  P-Ps  delay  from  near-field  data: 
2. 0-2. 5  seconds. 


♦Factor  multiplying  the  upgoing  waves  from  the  source. 
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Figure  46.  Comparison  of  synthetic  (heavy  lines)  and  observed  seismograms  at  HMNE  and  BFAK 
(Bluff,  Alaska). 


TABLE  12 


COMPARISON  OF  SYNTHETIC  ANO  08SERVED 
FOR  THREE  PAHUTE  MESA  E' 

h 

o_WAVE 

VENTS 

AMPLITUDES 

r* 

EVEMT/STATION 

b 

(nm) 

(s£c) 

mb 

r 

(nm) 

Tr 

(sec) 

mb 

MAST/HNME 

Observed 

594 

1.03 

5.99 

1159 

1.12 

6.25 

Synthetic 

Ratio 

524 

0.9 

0.92 

5.98 

906 

0.8 

1.07 

6.16 

FONTINA/HNME 

Observed 

819 

0.94 

6.17 

2014 

1.26 

6.43 

Synthetic 

Ratio 

1000 

1.2 

0.99 

6.25 

2397 

1.2 

1.28 

6.50 

CAMEMBERT/HNME 

Observed 

845 

0.99 

6.16 

1506 

1.09 

6.37 

Synthetic 

Ratio 

825 

1.0 

0.98 

6.16 

1453 

1.0 

1.15 

6.33 

MAST/BFAK 

Observed 

254 

0.99 

5.81 

466 

1.04 

6.05 

Synthetic 

Ratio 

321 

1.3 

0.97 

5.92 

656 

1.4 

1.04 

6.20 

FONTINA/BFAK 

Observed 

729 

0.99 

6.27 

1062 

0.95 

6.45 

Synthetic 

Ratio 

719 

1.0 

1.03 

6.24 

1532 

1.4 

1.20 

6.51 

camembert/bfak 

Observed 

375 

1.01 

5.97 

514 

1.06 

5.15 

Synthetic 

Ratio 

567 

1.5 

1.02 

6.14 

1040 

1.7 

1.09 

5.33 
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Figure  47.  Synthetic  seismograms  are  compared  to  observations  of 
PILEDRIVER  at  stations  HNME  and  BFAK.  The  observed 
records  have  been  hand  traced  from  film  chip  reproduc 
tions  (BFAK)  or  the  record  provided  with  the  SDL 
report  on  PILEDRIVER. 
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TABLE  13 


COMPARISON  OF  SYNTHETIC  ANO  OBSERVED 

FOR  PILEDRIVER 

P-WAVE 

AMPLITUDES 

STATION 

BFAK 

b 

(nm) 

(s£c) 

4 

c 

(nm) 

(sic) 

Observed 

65 

0.8 

5.31 

111 

0.8 

5.54 

Synthetic 

112 

0.78 

5.55 

153 

0.79 

5.69 

Ratio 

1.7 

1.4 

HNME 

Observed 

103 

1.0 

5.26 

210 

1.1 

5.53 

Synthetic 

184 

0.71 

5.66 

186  - 

0.70 

5.67 

Ratio 

1.8 

0.9 
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content  agreement  is  excellent  and  the  synthetic  seismogram 
amplitudes  are  about  40-70  percent  too  large.  In  other  words,  for 
PILEDRIVER  the  match  to  the  3FAK  data  is  about  the  same  as  for  the 
Pahute  Mesa  events. 

The  waveform  and  amplitude  comparison  at  HNME  is  rather  poor. 
Obviously,  the  relationship  between  the  HNME  and  BFAK  observations 
is  different  for  PILEDRIVER  than  it  was  for  the  Pahute  Mesa  events. 
This  can  be  seen  by  comparing  waveforms.  The  two  stations  are 
compared  in  terms  of  their  amplitude  ratios  and  their  magnitude  and 
period  differences  for  the  four  events  in  Table  14.  PILEDRIVER  is 
uniaue  in  having  the  HNME  record  be  much  longer  period  than  that  at 
BFAK.  If  the  instrument  response  were  indeed  as  we  think  it  was 
(Figure  45),  it  cannot  account  for  such  a  difference.  In  fact,  we 
would  expect  the  HNME  recording  (LRSM  seismometer)  to  be  shorter 
period  than  that  at  BFAK. 

What  can  we  conclude  from  the  comparison  at  these  two 
stations?  Our  PILEDRIVER  calculation  gives  synthetics  that  match 
the  3FAK  data  almost  as  well  as  did  our  synthetics  for  the  Pahute 
Mesa  events  which,  in  turn,  matched  the  HNME  data  even  better. 
However,  the  PILEDRIVER  waveform  comparison  at  HNME  is  not  good. 
There  seem  to  be  some  azimuthal  effects  not  present  in  the  Pahute 
Mesa  recordings. 

To  better  resolve  the  agreement  between  the  synthetic  and 
observed  short  period  P  waves  for  PILEDRIVER,  we  now  look  at  a  more 
extensive  data  set.  The  WWSSN  observations  of  PILEDRIVER  were 
collected  by  Hadley  and  wart  (1979).  These  data  are  displayed  in 
Figure  43,  which  is  reproduced  from  the  Hadley  and  Hart  report. 
Also  shown  in  the  figure  are  several  synthetic  seismograms  comouted 
from  our  two-dimensional  PILEDRIVER  source,  but  let  us  defer 
discussion  of  these  synthetic  records  until  we  have  described  the 
main  features  of  the  data. 

Hadley  and  Hart  (1979)  compared  the  observations  of  PILEDRIVER 
to  those  from  the  Pahute  Mesa  explosion  JORUM  at  ten  common  WWSSN 
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TABLE  14 


COMPARISON  OF  OBSERVED  DATA  AT  HNME  and  3FAK 

AMPLITUDES  MAST  CAMEMBERT  FONTINA  PILEDR  TVER 

b  (HNME/8FAK)  2.3  2.3  1.1  1.6 

C  (HNME/BFAK)  2.5  2.5  1.9  1.9 

MAGNITUDES 

mg  (HNME-8FAK)  0.18  0.19  -0.10  -0.05 

mg(HNME-BFAK)  0.20  0.21  0.31  -0.01 

PERIODS 

Tb  (HNME-BFAK)  0.04  -0.02  -0.05  0.2 

Tc  (HNME-8FAK)  0.08  0.03  0.31  0.3 
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stations.  The  amditude  ratios  are  listed  in  Table  15.  In  the 
northeast  quadrant  the  PILEDRIVER  amplitudes  are  markedly  smaller 
than  the  JORUM  amplitudes.  Hadley  and  Hart  (1979)  also  point  out 
that  at  the  northeast  stations  the  waveforms  for  P ILEDR I VER  are  much 
more  complex  than  those  for  JORUM.  Certainly,  we  can  see  in  Figure 
48  that  the  waveforms  from  the  stations  in  the  northeast  quadrant 
are  more  complex  than  most  of  the  others.  Hadley  and  Hart  conclude 
from  these  data  that  the  structure  beneath  P ILEDR I VER  is  introducing 
strong  azimuthal  effects. 

The  azimuth  from  P ILEDRI VER  to  HNME  is  60*,  within  the  region 
where  the  WWSSN  observations  are  relatively  complex.  However,  the 
HNME  seismogram.  Figure  47,  does  not  seem  to  be  specially  complex, 
so  the  effect  of  the  near  source  structure  (assuming  it  exists)  is 
not  as  strong  as  it  is  at  nearby  WWSSN  stations  like  WES,  OGD  or 
SCP.  We  did  see  (Table  14)  that  the  amplitude  of  P ILEDR I VER  at  HNME 
is  relatively  small  compared  to  the  3FAK  amplitude,  based  on  our 
experience  with  Pahute  Mesa  events,  though  the  amplitude  anomaly 
does  not  seem  to  be  as  large  as  it  is  for  the  WWSSN  stations.  All 
things  considered,  the  evidence  of  azimuthal  effects  is  enough  to 
reduce  the  usefulness  of  HNME  as  a  standard  for  comparison  of 
synthetic  and  observed  waveforms. 

We  now  compare  our  synthetic  seismograms  to  the  WWSSN  obser¬ 
vations  of  PILEDRIVER.  The  comparison  is  shown  in  Figure  ag  and  is 
done  taking  the  observations  as  a  whole,  rather  than  on  a  station- 
by-station  basis.  The  simplest  seismograms  were  written  at  KON, 
:<TG,  KIP,  NNA,  ARE,  MAT  and  LOR.  The  waveforms  at  those  stations 
are  very  similar  and  serve  as  a  standard  for  comparison.  They  are 
featured  by  a  b-phase  amplitude  that  is  60  percent  as  large  as  the  c 
onase.*  The  second  downswing  is  about  half  as  large  as  the  first, 
though  this  is  not  as  consistent  as  the  b/c  ratio.  Deviations  from 
this  basic  waveform  at  the  other  stations  are  nearly  always 
associated  with  lower  signal /noise  or  much  larger  coda  amplitudes, 

*  The  mean  b/c  for  these  seven  stations  is  0.61,  with  a  standard 
deviation  of  0.06. 
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TABLE  15 

AMPLITUDE  RATIOS,  PILEORIVER/ JORUM* 


Station  Azimuth  Ratio 

_  4  _ 


ARE 

133.1 

0.35 

KIP 

258.7 

0.50 

SHK 

309.1 

0.13 

COL 

336.1 

0.32 

AKU 

28.3 

0.25 

STU 

32.9 

0.27 

TOL 

46.0 

0.13 

WES 

67.0 

0.07 

OGD 

70.1 

0.06 

SCP 

71.2 

0.08 

133  <  4  <  33: 

Mean  .  0.30, 

standard 

deviation  »  0.12 

45  <  4  <  72: 

Mean  .  0.09, 

standard 

deviation  «  0.03 

AmDiituces  from  Figure  3.5.,  Hadley  and  Hart  (1979). 
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which  indicate  the  influence  of  structural  complexity.  An  important 
exception  is  COL,  which  happens  to  be  fairly  close  to  3FAK  where  the 
waveform  (Figure  47)  does  have  the  standard  pattern. 

The  synthetic  seismograms  were  computed  with  the  WWSSN  seis¬ 
mometer  response,  but  otherwise  are  as  in  Figure  47.  Seismograms 
are  plotted  in  Figure  48  for  three  values  of  t*,  which  are  0.8,  1.05 
and  1.3.  The  first  two  values  are  preferred  for  HNME  and  BFAK,  as 
explained  at  the  beginning  of  this  section.  The  largest  value,  1.3, 
is  the  average  t*  obtained  by  Helmberger  and  Hadley  (1979)  from 
simultaneous  fitting  of  near-  and  far-field  recordings  of  the  Pahute 
Mesa  events  JORUM  and  HANOLEY. 

The  synthetic  seismograms  are  not  in  good  agreement  with  the 
observed  waveforms  for  any  of  the  t*  values.  The  b/c  amplitude 
ratio  is  much  too  small.  Also,  the  secondary  arrival  that  we 
identified  with  spall  closure  in  the  calculation  (Section  5.7)  does 
not  seem  to  appear  on  the  observed  records.  Both  of  these 
deviations  from  the  observed  waveforms  are  consistent  with  there 
being  too  much  spallation  in  the  PIIEDRIVER  source  calculation.  We 
had  previously  remarked  that  this  was  the  case  when  examining  the 
near-field  velocity  data  (Section  4.1)  and  the  extent  of  cracking 
and  spallation  (Section  5.7). 

For  comparing  synthetic  and  observed  amplitudes,  we  essen¬ 
tially  follow  the  procedure  used  by  Hadley  and  Hart  (1979).  They 
list  the  maximum  peak-to-peak  amplitude,  corrected  only  for  the 
instrument  gain  at  1  Hertz.  They  then  corrected  the  data  to  a  range 
of  30’,  using  the  effective  geometric  soreading  factors  olotted  by 
Langston  and  Helmberger  (1975)  for  the  Jeffreys-3ul len  earth  model. 

In  Table  16  we  list  the  corrected  amplitude  data  from  Hadley 
and  Hart  (1979)  for  thirteen  selected  stations.  We  have  not 
included  the  North  American  stations  in  the  northeast  auadrant  where 
the  amplitudes  are  anomalously  low.  We  also  list  the  analogous  data 
for  the  b  phase  and  the  m^  computed  from  the  c  amplitude.  At  the 
bottom  of  Table  16  is  the  mean  m^  and  the  (logarithmic)  mean  b  and 
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TABLE  16 


AMPLITUDE  DATA  FOR  PILEDRIVER 


Station 

Range 

(Degrees) 

% 

Corrected 
c  Amplitude 

Corrected 
b  .amplitude 

ATL 

26 

5.63 

160 

110 

COL 

33 

5.94 

365 

146 

KIP 

40 

5.84 

648 

395 

KTG 

57 

5.46 

131 

89 

AKU 

60 

5.46 

137 

69 

NNA 

62 

5.90 

261 

144 

ARE 

68 

6.35 

748 

464 

KON 

73 

5.86 

329 

207 

LOR 

80 

5.65 

345 

173 

TOL 

81 

5.49 

180 

120 

STU 

82 

5.41 

135 

103 

SHK 

82 

5.76 

234 

164 

NAT 

86 

5.87 

384 

190 

Mean  m^:  5. 

74 

Mean  m 

i.  *  one  standard 
b 

deviation  = 

.  5.48  to  6.00 

Logarithmic 

Mean  c:  266 

Mean  c 

±  one  standard 

deviation  s 

149  to  475. 

Logarithmic 

Mean  b:  157 

Mean  b 

±  one  standard 

deviation  « 

91  to  270. 
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c  amplitudes.  Note  that  the  mean  m^  is  significantly  larger  than 
the  network  m,Q  values  we  mentioned  at  the  beginning  of  this  sec¬ 
tion.  These  other  values  may  be  biased  low  by  including  stations  in 
the  northeast  auadrant  with  anomalously  low  (compared  to  °ahute  Mesa 
events)  amplitudes. 

For  the  synthetic  seismogram  amplitudes  we  need  to  use  a  good 
"average"  value  for  the  geometric  spreading.  Since  most  of  the  data 
are  from  stations  beyond  40*  where  most  earth  models  give  about  the 
same  effective  1/R,  the  most  reasonable  procedure  is  to  use  the 
value  at  30*  for  the  Jeffreys-Bul len  model  used  to  correct  the  data 

r  w 

to  this  range.  This  value  is  8.4  x  10“°  km"1  (Langston  and 
Helmberger,  1975). 

The  amplitudes  for  the  three  synthetic  seismograms  plotted  in' 


Figure  48  are  listed  below 
dependent  instrument  response). 

(they  are  not 

corrected 

for  period 

t* 

b  Phase 

Tb 

Phase 

Tc 

0.8 

375 

0.76 

425 

0.78 

1.05 

167 

0.86 

163 

0.90 

1.3 

85 

0.97 

90 

1.39 

Since  the 

waveform  agreement 

is  not  good. 

the  most 

meaningful 

comparison 

is  for  the  b  phase. 

In  Table  16 

we  see  that 

this  mean 

value  is  157,  with  values  of  91  to  270  representing  one  standard 
deviation  from  the  mean.  Thus,  the  synthetic  with  t*  *  1.05  gives 
good  agreement  with  the  observed  amolitudes.  For  t*  =  0.3,  the 
amplitude  is  too  large,  and  it  is  too  small  for  t*  *  1.3.  As  far  as 
the  frequency  content  of  the  waveform  is  concerned,  the  t*  »  1.05 
case  seems  in  reasonable  agreement  with  the  data.  We  believe  the 
best  average  t*  for  NTS  explosions  recorded  at  teleseismic  stations 
is  near  1.0. 

Our  conclusion  is  that  the  P ILEDR I VER  source  has  about  the 
right  direct  P  amplitude  (the  b  amplitude  is  mainly  controlled  by 
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direct  P).  This  conclusion  cannot  be  quantified  with  mucn  precision 
because  of  uncertainties  about  the  best  average  values  for  the 
geometric  spreading  and  t*.  If  the  ? ILEDR I VER  source  is  in  error, 
our  general  impression,  based  on  many  poorly  defined  factors,  is 
that  it  is  too  large. 

The  two-dimensional  source  calculation  appears  to  include  too 
much  non-linear  interaction  with  the  free  surface.  The  later 
portion  of  the  synthetic  P  waveform  does  not  match  the  data, 
apparently  because  pP  is  too  greatly  suppressed  and  because  the 
spall  closure  energy  is  too  large,  or  is  timed  incorrectly. 

6.4  CONCLUSIONS 

A  comparison  of  this  kind  is  always  somewhat  subjective  and 
subject  to  qualification.  For  the  surface  waves  the  problem  is  how 
to  account  for  the  non-axi symmetric  component  that  is  present  in  the 
source.  If  we  believe  the  Rivers  and  von  Seggern  (1379)  solution 
for  this  component,  our  2-D  source  gives  surface  wave  amplitudes 
that  match  the  data.  Other  assumptions  about  this  component  of  the 
source  will,  of  course,  change  one's  perception  of  now  well  the 
surface  wave  data  are  matched.  Also,  the  body  wave  analysis  indi¬ 
cates  that  the  free  surface  interaction  effects  are  too  strong  in 
the  calculation,  at  least  at  body  wave  periods.  Reducing  the  free 
surface  interaction  to  better  match  the  P ILEDRI VER  snort-perioc  data 
would  probably  reduce  the  surface  wave  amplitude,  though  we  cannot 
estimate  the  amount  of  the  expected  reduction. 

for  body  waves  we  listed  our  conclusions  at  the  enG  of  tne 
last  section.  Our  best  estimate  is  that  tne  P  wave  source  ampliouce 
is  nearly  correct.  The  two-dimensional  calculation  ooes  not  apoear 
to  properly  model  the  free  surface  interaction  at  tne  short  perioas 
seen  on  the  P  wave  recordings,  apparently  because  free  surface 
interaction  effects  are  "overpredicted."  This  is  consistent  with 
the  comparison  of  computed  and  observed  near-fielo  data  (Section 
4.1),  which  indicate  that  the  calculation  included  too  much 
spallation. 
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VII.  COMPARATIVE  ANALYSIS  OF  ATI  AND  S-CUBED  GRANITE 

CALCULATIONS 


7.1  INTRODUCTION 

We  have  described  eleven  two-aimensional  calculations  of 

explosions  in  granite,  seven  aone  by  ATI  and  four  by  3-Cubed. 

Synthetic  body  and  su. face  wave  seismograms  have  been  computed  for 

all  these  calculations  and  we  have  listed  the  m.  ana  M  values. 

b  s 

In  this  section  we  present  a  unified  summary  of  these  results  so 
they  can  be  easily  compared. 

7.2  SUMMARY  OF  SURFACE  WAVE  MAGNITUDES 

In  Section  3.5  we  presented  the  Ms  for  the  ATI  granite 

calculations.  The  analogous  values  for  the  S-Cubea  calculations 

were  given  in  Section  5.5.  We  now  compare  the  two. 

The  same  path  model  (Table  2)  was  used  for  all  the  surface 

wave  calculations.  In  Section  6.2  we  pointed  out  that  this  path 

model  gives  M  values  that  are  about  0.35  higher  than  we  would 

expect  for  an  "average"  path  model.  If  we  are  to  compare  to  data, 

we  should  correct  for  tnis  overestimation. 

The  Mg  from  the  ATI  and  S-Cubed  granite  calculations  can  be 
directly  compared.  However,  the  properties  of  the  granite  in  the 
source  region  are  not  the  same  for  the  two.  Differences  in  elastic 
properties  lead  to  quite  predictable  differences  in  surface  wave 
propagation.  Bache,  Rodi  and  Harkrider  (1978)  point  out  that  for 
explosions  in  similar  materials,  the  surface  wave  amplitude  scales 
with  y '?  ,  where  u  is  the  shear  modulus.  Tne  ?  is  a  measure  or'  the 

30  30 

static  displacement  in  the  material. 

The  shear  modulus  for  the  ATI  granite  is  172  kbar,  compared  to 
u  »  206  kbar  for  the  S-Cubed  granite.  Thus,  if  all  other  factors 
were  equal,  we  would  expect  the  Ms  for  the  S-Cubed  calculations  to 
be  0.08  units  higher  than  for  the  ATI  calculations.  Differences 
other  than  this  may  be  attributed  to  differing  long  period  source 
levels. 
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A  good  estimate  for  the  source  amplitude  for  long  period 

surface  waves  is  provided  by  the  ifj  plotted  in  Figures  10  ana  22. 

The  amplitudes  are  all  scaled  to  a  common  yield,  0.02  Kt.  The 

W  |  3 

f?  I  at  15  seconds  for  the  ATI  calculations  is  about  4  m  for 
'  3 

253  m  depth,  about  8  m  for  531  m  depth,  ana  between  the  two  for 
the  others.  For  the  deep  (1  km)  S-Cubed  calculations,  the 
comparable  values  are  3.5  -  4.0  m^,  while  the  values  for  the 

3 

shallow  calculations  are  11  to  13  m  . 

The  for  all  the  granite  calculations  are  plotted  versus 
source  depth  in  Figure  49.  These  values  have  been  adjusted  to  an 
"average"  crust  oy  subtracting  0.35  from  the  M  presented 
earlier.  For  this  comparison  all  M  ,  including  the  observed  value 
for  PILEDRIVER,  have  been  scaled  to  150  kt  by  adding  log  (15Q/W). 

i  A  I 

As  expected  from  the  14'  I,  the  ATI  values  fall  between  the 

M  values  for  shallow  and  deep  S-Cubed  calculations. 

The  Ms  for  the  S-Cubed  calculations  show  a  remarkably 

consistent  dependence  on  depth.  As  shown  in  the  figure,  the 

-1  46 

parametric  dependence  is  about  h  *  .  This  is  a  depth  dependence 

that  is  very  much  stronger  than  expected  from  one-dimensional 

calculations  like  those  in  Figure  29.  It  may  be  compared  to  the 

0  33 

depth  dependence  of  h  J  in  the  semi-empirical  Mueller/Murphy 
model  (Section  3.4). 

In  Figure  50  we  plot  the  values  versus  yield.  The 
comparable  values  computed  with  the  S-Cubed  RDP  sources  469  ana 
472  (Figure  29)  are  also  shown.  This  is  simply  another  display  of 
results  that  nave  been  amply  discussed.  The  PILEDRIVER  observed 
value  is  quite  close  to  the  M  for  the  RDP  source  at  tne  appro¬ 
priate  depth,  but  is  0.32  units  lower  than  tne  for  tne  2-Q 
PILEDRIVER  calculation.  However,  recall  (Section  6.2)  that  a  study 
of  surface  wave  radiation  patterns  for  PILEDRIVER  by  Rivers  ana  von 
Seggern  (1979)  concluded  that  the  total  source  includes  a 
double-couple  contribution  (not  accounted  for  in  our  solution)  wnich 
acts  to  reduce  the  observed  PILEDRIVER  Rayleigh  wave  amplitudes  by  a 
factor  of  two  or  three.  The  double-couple  contribution  snoulo  be 
removed  for  consistent  comparison  with  the  theoretical  values  in 
Figure  50. 
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MAGNITUDE 


7.3  SUMMARY  OF  BODY  WAVE  MAGNITUDES 


Body  wave  magnitudes  were  computed  with  the  same  path  model 
for  both  the  ATI  (Section  3.6)  and  S-Cubed  (Section  5.6)  granite 
calculations.  In  comparing  these  values,  it  is  useful  to 
separate  effects  due  to  differing  elastic  properties  from  those  due 
to  source  coupling.  Bache,  et  al^  (1975)  point  out  that  the  P  wave 
amplitude  for  an  RDP  source  in  a  realistic  earth  model  is  roughly 

| Ap !  i  A  O  I 

proportional  to  a  \y  j  where  a  is  P  wave  velocity  and  Y  is  tne 
spectral  amplitude  of  the  reduced  velocity  potential  at  the 

appropriate  period.  We  computed  an  effective  RDP,  ’V  ,  for  the 

two-dimensional  S-Cubed  calculations  in  Section  5.9  ano  found  that 
it  was  not  too  different  from  the  RDP  from  a  one-dimensional 
calculation  in  the  same  material,  especially  for  the  deeper  events. 
The  same  is  probably  true  for  the  ATI  calculations. 

The  P  wave  velocities  for  the  ATI  ana  S-Cubed  calculations  are 
4.403  km/sec  and  5.35  km/sec.  Thus,  if  the  effective  source 
coupling,  Yg  ,  were  the  same  for  the  two,  we  would  expect  the  mb 
for  the  S-Cubed  calculations  to  be  0.08  units  higher.  Note  that 
this  is  the  same  as  the  expected  M$  difference  due  to  elastic 

properties. 

We  will  be  plotting  the  ATI  and  S-Cubed  body  wave  magnitudes 
versus  source  depth  and  yield.  We  want  to  include  the  observed  data 

on  the  same  plots,  so  should  be  sure  that  the  theoretical  values 

include  our  best  estimate  of  the  path  parameters.  In  Section  6.3  we 
compared  observed  and  computed  values  for  PIIEDRIVER  as  well  as  we 
could  and  concluded  that  our  two-dimensional  calculation  gives 
seismograms  that  are  probably  about  the  right  amplitude. 

We  will  now  plot  the  ATI  ano  S-Cubea  body  wave  magnitudes 

versus  source  depth  and  yield.  We  also  wish  to  include  the  observed 
data  on  the  same  plots.  The  synthetic  seismograms  presented  in 

Sections  III  and  V  did  not  include  our  "best"  estimate  for  the 

average  attenuation  and  geometric  spreading,  so  these  mb  data  are 

shifted  (by  a  constant)  from  the  values  we  consider  suitable  for 

comparison  to  the  data. 
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What  about  the  data?  Alewine,  Young,  Springer  and  Klepinger 
(unpublished  report)  have  given  m^  values  for  the  granite  ex¬ 

plosions  PILEDRIVER,  SHOAL  and  HARDHAT.  These  values  are  comparable 
to  our  mF.  Based  on  the  ratio  of  b  to  c  amplitudes  for  PILE- 

DRIVER  (Section  6.3),  we  estimate  that  a  b  phase  magnitude  (mF) 
would  be  about  0.25  units  smaller. 

The  Alewine,  et_  a]_.  m£  for  PILEDRIVER  is  5.47.  In  Table 

16  we  saw  that  the  mean  for  thirteen  WWSSN  stations  was 

5.74.  The  b  amplitude  for  our  "best"  synthetic  PILEDRIVER  seismo¬ 
gram  was  in  good  agreement  with  the  b  amplitude  from  these 
stations.  We  do  not  believe  the  lower  magnitude  of  Alewine,  et  al . 
gives  a  more  accurate  estimate  for  the  teleseismic  p  wave 
amplitudes,  so  we  feel  justified  in  normalizing  all  the  magnitudes 
by  a  factor  that  causes  the  synthetic  and  observed  m^  for 
PILEDRIVER  to  plot  close  together.  For  plotting  observed  values  for 
all  three  granite  explosions,  it  is  more  convenient  to  use  the 
Alewine,  et  al_.  magnitudes,  so  we  obtain  a  PILEDRIVER  m^  by 

subtracting  0.25  from  the  Alewine,  et  al.  mF  of  5.47.  The 

result  is  5.22.  If  we  subtract  0.3  from  the  synthetic  for 
PILEDRIVER  (Table  10),  it  becomes  5.32,  and  the  0.1  difference  is 
about  right  for  indicating  the  agreement  between  theoretical  and 

observed  values.  Therefore,  we  plot  the  Alewine,  et_  ajL  observed 
magnitudes  and  al  1  theoretical  magnitudes  are  adjusted  by  sub¬ 
tracting  0.3  units  before  plotting.  Alternatively,  we  could  have 
added  0.3  to  the  Alewine,  et_  al_.  observed  magnitudes  to  make  them 
(at  least  PILEDRIVER)  consistent  with  the  observed  WWSSN  magnitudes. 

In  Figure  51  we  plot  the  adjusted  oody  wave  magnitudes  for  the 
two-dimensional  calculations  versus  source  depth.  The  S-Cubed 
calculations  and  the  observed  m^  have  been  scaled  to  150  kt  by 

adding  log  (150/W).  That  is,  we  scale  without  accounting  for  any 
yield-dependent  freauency  shift.  The  dependence  on  depth  of  the  "b" 
phase  magnitude,  m^,  is  auite  weak  for  both  the  ATI  and  S-Cubed 
calculations.  The  "c"  phase  magnitude  is  more  sensitive  to  free 
surface  effects,  but  it  too  has  no  strong  depth  trend. 
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Figure  51. 


The  values  for  all  the  two-dimensional  calculations  are 
plotted  versus  source  depth.  The  S-Cubed  calculations  at 
the  same  yield  are  connected  with  a  dashed  line.  The  ob¬ 
served  values  for  HARDHAT,  SHOAL  and  PILEDRI VER  are  also 
plotted  with  m,  =  m?  -  0.25.  All  magnitudes  have  been 
scaled  to  150  KT.  0 
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The  S-Cuoea  calculations  show  much  stronger  coupling  into 
short  perioc  P  waves  than  the  ATI  calculations.  Accounting  for  tne 
elastic  propagation  effect,  we  estimate  that  the  equivalent  RDP  for 
the  S-Cubed  calculations  is  about  2.2  times  larger  near  1  Hertz  than 
that  for  the  ATI  granite.  Since  tne  long  period  level  of  the 
equivalent  ROP  for  the  ATI  calculations  is  between  the  values  for 
the  shallow  and  deep  S-Cubed  calculations,  this  means  that  the 
effective  ROP  source  function  for  the  S-Cubed  granite  is  much  more 
peaked. 

The  body  wave  magnitudes  are  plotted  versus  explosion  yiela  in 
Figures  52  and  53.  As  pointed  out  in  Section  V,  the  is 
about  the  same  for  the  one-  and  two-dimensional  S-Cubed  source 
calculations.  The  two-dimensional  effects  are  apparent  in  the 

m^,  especially  for  the  shallow  events. 

Our  final  plot,  Figure  54,  shows  the  residual  as  a 

function  of  source  aepth.  We  also  plot  this  residual  for  the 

Muel ler/Murphy  RDP  source.  The  difference  between  the  S-Cubed  ana 
ATI  calculations  is,  once  again,  tne  strong  dependence  of  on 
source  aepth  for  the  S-Cubed  calculations.  The  residuals  for  the 

Muel ler/Murphy  RDP  source  have  a  similar  depth  dependence,  but  for 
quite  different  reasons.  For  an  RDP  source  in  an  elastic  earth 
model,  the  m^  is  strongly  dependent  on  depth  because  of  the 

interference  of  P  and  pP. 
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Figure  5 
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The  b  phase  magnitude,  mS,  is  plotted  versus  yield. 


The 


b> 

values  for  the  ATI  calculations  fall  within  the  indicated 
range.  The  observed  values  are  the  published  mP  minus 
0.25  for  HARDHAT,  SHOAL  and  PILEDRIVER.  0 
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Figure  52.  The  c  phase  magnitude,  mjj,  is  plotted  versus  yield. 

lished  (network  average)  magnitude  values  for  HAROHAT 
SHOAL  and  PIISDRIYER  are  also  slotted. 
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The  difference  between  body  and  surface  wave  magnitude  is  plotted 
versus  source  depth. 
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APPENOIX  A 


!  SYNTHETIC  SEISMOGRAMS 

FROM  COMPLEX  SOURCE  CALCULATIONS 

A. 1  INTRODUCTION 

There  is  increasing  interest  in  the  use  of  large  scale  finite 
difference  or  finite  element  computer  programs  to  simulate  the 
complex  nonlinear  material  behavior  that  occurs  in  the  region 
immediately  surrounding  earthquake  or  explosion  sources.  To  fully 
appreciate  the  significance  of  the  output  of  these  programs,  it  is 
important  to  be  able  to  compute  the  ground  motions  at  ranges  much 
larger  than  the  maximum  dimension  of  the  computational  grid.  As 
long  as  the  material  response  can  be  assumed  to  be  linearly  elastic 
outside  the  immediate  source  region,  analytical  techniques  can  be 
used  to  continue  the  wave  field  to  the  ranges  of  interest.  In  this 
Appendix  we  describe  these  techniques. 

When  the  source  is  assumed  to  be  embedded  in  an  elastic 
wholespace,  an  "equivalent  elastic  source"  representation  can  be 
constructed  by  expanding  the  outgoing  wave  field  in  spherical 
harmonics  (e.g.,  Bache  and  Harkrider,  1976).  The  source  is  then 
expressed  in  terms  of  multipolar  coefficients  that  can  be  used  as 
input  to  programs  for  computing  synthetic  seismograms.  This 
technique  has  been  used  to  study  the  results  from  complex  explosion 
(Cherry,  et.  al_. ,  1975,  1976;  Bache  and  Masso,  1973)  and  earthquake 
(Bache,  et.  al . ,  1980)  sources  computed  with  finite  difference 
programs. 

When  boundaries  (e.g.,  a  free  surfac.)  are  present  in  the 
vicinity  of  the  source,  the  ground  motions  at  ranges  outside  the 
numerical  grid  can  be  computed  with  a  wave  field  continuation 
method.  Based  on  an  elastodynamic  representation  theorem  (e.g., 
Burridge  and  Knopoff,  1964).  This  method  requires  the  convolution 
of  stresses  and  displacements  from  the  elastic  region  of  the 
numerical  source  calculation  with  Green's  functions  representing 
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the  elastic  wave  propagation  through  the  medium  outside  the  source 
region. 

In  subsequent  sections  we  begin  with  the  formulation  of  the 
basic  theory.  For  earthquake  sources  the  wave  field  continuation 
must  be  done  in  three  dimensions.  However,  many  important  complex 
explosion  phenomena  are  axisymmetric  and  the  theory  is  specialized 
for  this  case.  Some  of  the  development  for  the  axisymmetric  case 
can  also  be  used  for  three-dimensional  sources  in  plane-layered 
earth  models,  since  the  Green's  functions  are  azimuth-independent 
for  such  structures. 

The  details  of  the  implementation  depend  on  the  choice  of 
Green's  functions  for  the  propagation.  Our  emphasis  here  is  on  the 
Rayleigh  wave  modes  for  plane-layered  earth  models.  Test  problems 
which  illustrate  this  case  and  demonstrate  the  accuracy  that  can  be 
achieved  are  presented.  In  the  body  of  the  report.  Sections  III  and 
V,  we  discuss  synthetic  far-field  P  wave  seismograms.  For  these 
synthetic  seismograms  we  used  Green's  functions  computed  with  the 
propagator  matrix  method  of  Fuchs  (1966). 


A.  2  FORMULATION 


For  the  analytic  continuation  we  use  an  elastodynamic 
representation  theorem  (Burridge  and  Knopoff,  1964).  Assuming  no 
body  forces  are  present,  this  theorem  may  be  written 


“f 


G1.*^ 
J  J 


dA 


(A.l) 


where  the  repeated  indices  imply  summation  and  the  following 
notation  has  been  used: 

uF(x„,t)  :  Components  of  the  displacements  at 

Xp  outside  the  source  region; 
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S  .  ( X«t»  X-,  t  ) 

j  ~r 


S]  (xM,xF,t): 
J  k  * 


rM(x, 


t) 


M 
u  . 
J 


Displacement  in  the  direction  j  at 
xM  due  to  a  unit  impulsive  force 
at  xF  in  the  direction  i; 

The  components  of  the  stress  tensor 
at  x„  due  to  a  unit  impulsive 

n 

force  at  xc  in  direction  i; 

r*  f"  ** 

Components  of  the  monitored 
traction  vector  on  the  monitoring 
surface; 

Components  of  the  monitored 
displacement  vector  on  the 
monitoring  surface; 

Components  of  the  normal  to  the 
monitoring  surface,  with  the 
positive  sense  being  for  vector 
pointing  away  from  the  source; 


S 


M 


Area  of  the  monitoring  surface. 


Also,  convolution  is  indicated  by 


f*g 


f (t-r)g(T) 


dx 


(A.2) 


Thus,  the  solution  at  some  distant  point  is  obtained  by  first 
convolving  suitaole  Green's  functions  with  stresses  and 
displacements  monitored  on  some  surface  that  encloses  the  region  of 
inelastic  material  deformation.  The  final  solution  is  then  obtained 
by  integrating  the  convolution  products  over  the  monitoring 
surface.  The  restrictions  on  this  technique  are,  first,  that  the 
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region  outside  the  monitoring  surface  must  be  linearly  elastic. 
Second,  if  the  problem  domain  is  bounded,  the  Green's  functions  must 
satisfy  appropriate  boundary  data.  The  elastic  medium  may  be 
anisotropic  or  have  discontinuous  properties,  though  the  computation 
of  Green's  functions  is  much  more  difficult  for  such  media  than  for 
homogeneous  isotropic  media. 

To  this  point  the  formulation  is  entirely  general.  Let  us  now 
consider  cases  where  the  source  and  propagation  medium  are 

axisymmetric.  This  is  the  appropriate  geometry  for  explosions  near 
the  free  surface  (or  other  horizontal  boundaries)  in  plane-layered 
earth  models. 

As  shown  in  Figure  A.l,  stresses  and  displacements  are 

monitored  on.  some  surface  with  a  typical  point  M.  We  wish  to 

compute  the  displacement  field  at  F.  Let  the  monitored  quantities 

at  M  be  denoted  by  o*?„  aM„,  <?£!,,  uM  and  uM.  Then  the 

rr  rz  zz  r  z 

components  of  the  traction  vector  and  displacement,  referred  to  a 
local  coordinate  system  (x^,  x^,  x^)  at  F,  are  given  by: 

T?  »  (aj?r  cos©  +  a^z  sin©)  cos  S  , 

T2  *  '  ^arr  C0S&  +  arz  sine)  sin  ’  ’ 

_M  M  M 

t3  *  ap2  COS©  +  azz  Sine  , 

M  M 

'Jj  =>  >Jr  COS  $  ,  (A. 3) 

M  VI 

u,  =  -  ur  sin  ;  , 
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At  the  top  we  show  the  inelastic  region  and  monitoring 
surface  for  a  typical  axi symmetric  problem.  At  a  moni¬ 
toring  point  M  the  outward  normal  to  the  monitoring  sur¬ 
face  is  n.  The  coordinate  system  for  defining  the 
Green's  function  Is  shown  in  the  horizontal  plane  (viewed 
from  positive  X3)  at  the  bottom. 
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=  sing  . 


n,  =  coso  cos  S  ,n,  =  -cos&  sin  -  , 


Also,  the  displacements  at  F  in  the  global  R-z  system  are  related  to 
the  displacements  in  the  x^,  x^,  x^  system  by 


u^  =  u^  cos  u,  —  u ^  sin  y  , 


(A.A) 


«,F  -  A 


Using  these  coordinate  transformations,  (A.l)  may  be  written: 


r  (G 


„R*  M  ^  -R*  M 


R  4...M  ^  CR  *  M, 


rV  +  G‘  V  )  -  (S"  *u"  +  S*  *u  )  dz 
r  rr  z  r z‘  v  rr  r  rz  z' 


/  f/rR*  M  .  «R+  M  ,  /CR  .  eR  +  M,  . 

r  (Gz*o22  ♦  Gr*«r2)  -  (Sr2*ur  *  S22*u2)  dr  , 

JM  L 


-  r  (GZVM  +  GZ*<7M  ) 
/  '  r  °rr  z  rz' 


(Sz  *uM  +  Sz  *u  )  dz 
'  rr  r  rz  z ' 


/„  \,  -Z*  M  .  _z*  M  «  ,,-z  M  ^  M  .  <-z  +  M,  . 

r  [(Gz  azz  +  V^rz)  -  (Srz  ur  ^r  +  Szz  uz}  dr  • 

R  R 

where  the  quantities  Gr,  Gz,  etc.,  are  azimuthal  averages  of  the 
Green's  functions  in  the  x^,  x^,  x^  system.  The  explicit 

formulas  for  the  azimuthal ly  averaged  Green's  functions  are  the 
following: 


Gj  COS';/  COS  ; 


1  2 

-  G,  sin  (£+>)  +  sinking  cd 


1  2 

G^  cos g  -  G^  sing  dd  , 


(A. 6) 
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G3  cos 


-  G^  sin  d i  , 


G^  d  <t>  , 


(A. 7) 


where  G^(xM,Xp,t)  ,  etc.,  are  evaluated  at  x  =  0.  Also, 


iS,v|  m  (X^»X2»X3) ,  with 

X1  *  RF  +  rM  *  2rM  Rf  cos*‘ 


x2  -  0  , 


X3  “  ZM  “  2F  * 


(A.8) 


For  axisymmetric  problems  formulated  to  compute  the  solution  in  the 
«j>,  k  domain  (e.g.,  modal  solutions),  the  azimuthal  integrations  in 
(A. 6)  and  (A. 7)  may  be  computed  analytically.  In  the  p,  t  domain 
(e.g.,  generalized  ray  theory),  these  integrations  cannot,  in 
general,  be  performed  analytically.  However,  in  the  far-field,  an 
analytic  approximation  for  these  integrals  may  be  used. 

The  formulas  for  The  azimuthally  averaged  stress  Green's 

D 

functions,  Srr,  etc.,  are  given  by: 

2ir  (S^  COS 2£  -  sin2?  +  sin2?  cost 

3rr  a  /  at>’ 

^  ~  ''Sil  cos2’  *  S2i  sin2’  +  S22  sir|2^  )  sin- 


Zl  (S^  cos  £  -  sin,-)  cos- 


0  1  -  (S?,  cos^  -  S?9  sinC)  sint 


(A. 9) 


2  TT 

74 


S2n  cos 2C 


-  sh  sin2 ?  +  S L  sin2£ 
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S31  cos^-  S3 


CTt 

Szz  =  [  S33  cosc  ‘  S33 

J  0  L 


S33  d*‘ 


It'  the  azimuthal  averages  in  (A. 6)  and  (A. 7)  can  be 
analytically  determined,  it  is  not  necessary  to  use  (A. 9).  Instead, 
the  constitutive  equation  of  the  medium  at  the  monitoring  point  can 
be  used.  Thus,  for  an  isotropic  medium 


SA  -  x 
rr 


+  2»tf 


/  3GZ  GZ  3G2  \  3GZ 

Sz  ,  x  _£  +  _n  +  _i.  )+2u— z-  , 
°zz  \  ar  r  az  /  az  ’ 


\  az  ar  ', 


(A. 10) 


where  x  and  w  are  Lame's  constants,  and  formulas  of  the  same  type 

can  be  written  for  SR  .  and  . 

rr  zr  zz 

The  basic  formulation  for  analytical  continuation  of  axi- 
symmetric  source  computations  is  now  complete.  These  equations  are 
valid  for  a  monitoring  surface  of  arbitrary  shape.  Implementation 
requires  the  computation  of  Green's  functions  appropriate  for  the 
application  of  interest. 
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A. 3  FORMULATION  FOR  SURFACE  WAVES 

For  many  problems  the  primary  interest  is  in  the  normal  mode 
contribution  to  the  solution  for  waves  in  layered  media.  For 
example,  the  far-field  fundamental  mode  Rayleigh  wave  is  commonly 
observed  and  studied  for  seismic  events.  For  such  problems  the 
analysis  is  more  conveniently  done  in  the  frequency  -  wave  number 
domain. 

We  will  develop  the  solution  using  Green's  functions  computed 
following  Harkride*  (1964,  1970),  though  the  theory  is  easily 

modified  to  use  alternative  techniques  for  computing  Rayleigh  wave 
modes. 

For  the  azimuthal ly  averaged  Rayleigh  wave  displacement. 
Green's  functions  for  a  particular  mode  are: 

Gr  *  ilT  ~R  *0  HpVR)  . 

G2  *  ’*  A)  *0  “(z)  J0(kr)  Hl2)(kR>-  (A. 11) 

G l  -  *  X(2)  Jj(kr)  H<2)(kR)  , 

G(  *  '  4r  “(I)  Jo(kl')  Ho2)(kR)  • 

where  X(z)  and  W(z)  are  the  normalized  horizontal  ano  vertical 
displacement  eigenfunctions  and  eQ  is  the  ellipticity.  The  A^  is 
the  amplitude  response,  k  is  wave  number,  r  is  the  radial  distance 
to  a  point  on  the  monitoring  surface  and  R  is  the  range. 

The  stress  Green's  functions  are  easily  computed  from  (A. 10) 
using  the  following  identities  given  by  Harkrider  (1970): 
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where  t  (z)  and  Z(z)  are  the  normalized  horizontal  and  vertical 
stress  eigenfunctions. 


A. 4  TEST  CALCULATIONS 

To  illustrate  the  method  and  determine  its  accuracy,  we 
compare  exact  and  analytically  continued  Rayleigh  waves  for  a 
spherically  symmetric  explosion  source  which  is  representea  by  the 
potential  plotted  in  Figure  A. 2.  The  exact  solution  is  computed 
with  the  method  of  Harkrider  (1964).  The  monitored  solutions  for 
input  to  the  wave  field  continuation  method  were  computed  in  two 
ways,  analytically  with  generalized  ray  theory  and  numerically  with 
the  SIm IS  finite  element  program  (Frazier  and  Peterson,  1974).  The 
geometry  for  the  generalized  ray  theory  calculations  of  the  stresses 
and  displacements  at  ten  evenly  spaced  monitoring  stations  is  Shown 
in  Figure  A. 3.  The  elastic  properties  of  the  halfspace  are  a  *  6.0 

3 

km/sec,  e  *  3.5  km/sec,  p  *  2.7  gm/cm  .  However,  the  monitored 
quantities  can  be  scaled  to  other  materials  with  the  same  Poisson's 
ratio  by  dividing  the  displacements  by  the  shear  modulus,  »,  and  the 
time  by  the  velocity,  a.  This  scaling  is  easily  seen  from  the 
reduced  displacement  potential  representation  for  an  explosion.  The 
displacement  is  related  to  this  potential  by 


u 


Ijt  -  R/a) 

3 


(A. 13) 


where  R  is  the  spherical  coordinate.  For  a  pressure  profile  applied 
to  the  surface  of  a  spherical  cavicy  of  radius  R  the  Fourier 
transformed  reduced  displacement  potential  is  (Mueller  and  Murphy, 
1971) 
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-  7 

10  * 

FREQUENCY  (Hz) 


Figure  A. 2.  The  spectral  amplitude  of  the  Fourier  transform  of  the 
reduced  velocity  potential  used  for  the  explosion  simu 
lations. 
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SURF  AC 


Figure  A. 3.  The  geometry  for  the  generalized  ray  theory  calculation 
for  a  spherically  symmetric  explosion  in  a  halfspace. 
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(A.  14) 


.1,  ,  '3(“>  Re  «2 

=  - -  3  T-?  p  O 

y  4a  +  4iuaRg  -  a  id  R^  /S 

where  P(w)  is  the  transformed  pressure  time  history.  Fixing  P(u) 
and  Rg,  the  simple  scaling  is  obtained. 

In  Figure  A. 4  we  compare  the  surface  wave  computed  with  this 
method  with  the  exact  solution  for  two  propagation  media.  The  first 
is  the  elastic  halfspace  of  the  source  calculation  with  Q  =  50. 
The  second  comparison  is  for  a  realistic  earth  model.  Bache,  et_. 
al .,  (1973),  pointed  out  that  there  is  often  a  need  to  account  for 
the  material  properties  in  the  immediate  source  vicinity  being 
different  from  the  average  properties  along  the  travel  oath  and 
presented  an  approximate  techniques  for  accounting  for  this 
difference.  An  analogous  approach  was  used  for  the  second 
comparison  in  Figure  A. 4.  The  source  region  was  represented  by  the 
NTS/TUC  model  of  Bache,  et  aJL,  with  the  top  two  kilometers  replaced 
by  a  layer  with  a  =  3.5  km/sec,  s  =  1.75  km/sec,  p  =  2.0  am/cnT. 
The  monitored  solutions  were  scaled  to  be  appropriate  for  this 
material.  Then  the  source  excitation  terms  in  the  Green's  functions 
were  computed  with  this  modified  structure,  while  the  Hankel 
functions  representing  the  propagation  were  compu+  ■  with  the  phase 
velocities  for  the  original  NTS/TUC  model.  The  transition  between 
the  source  region  and  propagation  path  models  was  accounted  for  by 
the  transmission  coefficient 


i 


where  subscripts  1  and  2  represent  the  source  and  propagation 
models,  respectively. 


The  exact  seismograms  are  in  very  close  agreement  with  those 
from  the  wavefield  continuation  method.  In  Figure  A. 5  we  compare 
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NTS/TUC 


Figure  A. 4.  Comparison  of  exact  (top)  seismograms  with  those  com¬ 
puted  with  the  wavefleld  continuation  method.  The 
numbers  at  the  left  are  peak  amplitudes  in  microns. 

The  range  was  1,000  km  for  the  halfspace  and  728  km 
for  the  NTS/TUC  example.  The  seismograms  were  filtered 
by  the  WWSSN  15-100  instrument  response. 
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the  two  in  the  frequency  domain.  We  see  that  the  errors  are  less 
than  10%  over  the  entire  band  from  0.5  to  40  seconds.  After  scaling 
the  halfspace  comparison  to  the  source  material  of  the  NTS/TUC 
model,  the  ratio  for  the  two  cases  is  nearly  the  same  for  periods 
less  than  10  seconds.  At  longer  periods  the  contributions  from  the 
integrations  (A. 5.)  along  the  side  and  bottom  of  the  cylindrical 
monitoring  surface  have  opposite  sign  and  are  nearly  the  same  size. 
Then  small  errors  in  the  static  values  of  the  monitored  solutions 
(which  were  only  carried  out  to  1  second)  are  weighted  by  different 
Green's  functions  for  the  two  cases,  leading  to  different  errors. 
At  shorter  periods  the  primary  reasons  the  exact  and  wavefield 
continuation  solutions  diverge  are  inaccuracies  in  the  monitored 
solutions  and  the  lack  of  resolution  of  the  sparse  spatial  sampling. 

A  test  problem  more  closely  related  to  the  ultimate  problems 
of  interest  was  done  with  the  SW IS  finite  element  program  (Frazier 
and  Peterson,  1974)  using  the  geometry  shown  in  Figure  A. 6.  The 
grid  dimensions  were  50  by  50  meters.  An  exact  simulation  of  an 
explosion  source  requires  application  of  a  pressure  to  the  surface 
of  a  spherical  cavity.  However,  in  this  geometry  it  was  more 
convenient  to  apply  the  pressure  inside  a  cylindrical  cavity  with  50 
meter  radius  and  height.  7ne  stresses  and  displacements  were 
monitored  at  the  indicated  20  stations  and  were  convolved  with  the 
reduced  displacement  potential  of  Figure  A. 2.  The  amplitudes  were 
normalized  by  requiring  the  product  of  cavity  pressure  and  volume 
for  the  cylindrical  cavity  to  be  the  same-  as  for  a  spherical  cavity 
consistent  with  the  reduced  displacement  potential. 

The  seismograms  for  the  a  =  6.0  '<m/sec  halfspace  are  compared 
in  Figure  A. 6.  The  spectral  comparison  is  also  shown.  The  com¬ 
parison  is  done  for  a  solution  using  all  twenty  monitoring  stations 
and  for  one  with  alternate  modes  discarded.  With  20  stations,  the 
errors  are  less  than  10  percent  from  40  seconds  to  0.3  seconds,  but 
the  errors  are  much  larger  at  long  periods  with  fewer  stations. 
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Figure  A. 6,  Schematic  view  of  the  axisymmetric  grid  for  the  SWIS 
calculation.  The  twenty  monitoring  solutions  are  In¬ 
dicated. 
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Figure  A. 7.  Temporal  and  spectral  ratio  comparison  of  exact  and  analytically  continued  halfspace 

Rayleigh  waves  after  propagation  to  1000  km  and  filtering  by  the  WWSSN  15-100  instrument 
response. 


A. 5  ANALYSIS  OF  THE  LONG  PERIOD  SOLUTION 


The  test  problems  presented  in  the  previous  section 
demonstrate  that  the  technique  is  quite  accurate  over  most  of  the 
period  range  of  interest.  However,  at  long  periods  the  errors  can 
be  fairly  large,  even  for  these  ideal  test  problems  where  the  late 
time  values  of  the  monitored  solutions  are  defined  with  considerable 
accuracy. 

In  the  highly  nonlinear  complex  source  calculations  of 
practical  interest  it  is  usually  impossible  to  continue  until 
motions  entirely  cease  at  the  monitoring  surface.  These 
calculations  generally  must  be  terminated  after  a  few  seconds  when 
the  motions  are  small,  but  finite.  Therefore,  the  static  stresses 
and  displacements,  which  are  what  determine  the  solution  at  periods 
much  larger  than  the  duration  of  the  calculation,  will  not  be 
determined  with  great  precision. 

The  test  problems  were  done  on  a  sparse  spatial  grid  and  much 
of  the  long  period  error  may  be  attributed  to  inaccurate  spatial 
quadratures  (Figure  A. 7).  This  inaccuracy  can  be  nearly  eliminated 
by  monitoring  at  each  node  in  the  discrete  grid  of  the  source 
calculation.  A  more  serious  concern  is  with  the  effect  of  errors  in 
the  computed  static  solutions,  which  are  generally  the  values  of  the 
monitored  quantities  at  the  last  time  step. 

Let  us  examine  the  low  frequency  limit  of  the  vertical 
Rayleigh  wave  for  axisymmetric  problems.  Assume  that  the  monitoring 
surface  is  a  cylinder  of  radius  a  and  depth  b.  The  monitored 
stresses  and  displacements  are  assumed  to  reach  some  static  value. 
In  this  case  their  transforms 


aM  aM 

azz’  CTrz’ 


are  proportional  to  at  low  frequency.  Then  with  some  algebra 
it  can  be  shown  that 


•-« 

. 

; 


‘I 


185 


1 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


1 im  aF,  > 
»»0  U2(r) 


•f*&  «o(2)  <kr> 


f  »22  (r’l>)  r0  S 


•7 


+  /  a  djr(a,z)  dz 


-kc0 


26  2 

1" 
a 


aM  r  , 
°zz  r0  dr0 


I  co 


a 

/ 


aM  wv 

r0  azr  dr0 


(A. 16) 


+  k  e  a 
o 


/[ 


y(3x  +  2u)  aM  £  aM 
x  +  2U  r  ■  2  ffrr 


-  z 


U" 

J 

!  1 

1 +  °| 

2 

in 

/  zr 

) 

. 

where  all  temporal  quantities  have  been  Fourier  transformed  and 
depend  on  <u.  The  x,  y  are  the  Lamfe  constants. 

The  first  two  terms  in  the  brackets  are  simply  the  applied 

A 

vertical  force,  F2(u>).  That  is. 


A 

C 


'z 


(«) 


a  a. 


(a,z)  dz  .(A. 17) 


For  most  source  calculations  the  medium  is  initially  at  rest  and  o 
external  or  body  forces  are  applied  (the  effect  of  gravity  can  be 
subtracted  from  the  monitored  solutions).  If  so,  conservation  of 
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Figure  A. 8.  The  Rayleigh  wave  spectral  amplitude  for  the  halfspace 
example  of  Figure  A. 4  Is  plotted  together  with  the  con¬ 
tributions  from  the  monitored  displacements,  the  stresses 
associated  with  vertical  tractions  (C2),  and  the  stresses 
associated  with  radial  tractions  (Cj. 

P 
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linear  momentum  requires  that  the  vertical  force  ana  the  vertical 
impulse  vanish  at  static  equilibrium.  That  is, 

• 0  • 

■  111  /  *-«• 
yo 

Then  the  Fz(«)  terms  in  (A. 16)  are  proportional  to  u>  at  low 
frequency  while  all  other  terms  are  0(1). 

Also,  if  F  (t)  approaches  zero,  but  I  (t)  does  not,  then 

•  A  * 

the  terms  proportional  to  Fz(«)  in  (A. 16)  are  the  same  order  as 
the  others.  Further,  if  F  (t)  i  0  at  late  time,  this  contribution 
is  proportional  to  1 /«  and  dominates  the  low  frequency  solution. 

When  (A. 18)  is  satisfied,  the  low  frequency  limit  is 


(A. 19) 


Interpretation  of  this  result  is  more  convenient  when  it  is  written 

AM  PM 

in  terms  of  tractions,  Tr,  T  ,  instead  of  stresses  ana  for 
a  monitoring  surface  of  aroitrary  shape.  Then 

F>)  -  /  Tj  (r,z)  dA  ,  (A. 20) 

•4 
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and  (A. 19)  becomes 
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(A. 21) 


Thus,  the  low  frequency  behavior  of  the  solution  depends  only  on  the 
vertical  and  radial  tractions  and  the  radial  displacements. 

As  an  illustration  of  the  relative  importance  of  the  various 
terms  comprising  the  solution,  we  examine  the  halfspace  example  of 
Figure  A. 4.  The  Rayleigh  wave  spectral  amplitude  (before 
application  of  the  filters  for  the  0  and  seismometer)  is  plotted  in 
Figure  A. 8  together  with  the  spectra  for  three  contributions  that 
make  up  the  total  solution.  In  this  example  the  displacement 
contributions  dominate  at  long  periods  and  the  radial  traction  terms 
are  relatively  small  over  the  whole  frequency  band.  However,  the 
relative  importance  of  these  terms  is  dependent  on  details  of  the 
problem  at  hand. 
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APPENDIX  B 

AN  "EQUIVALENT"  POINT  SOURCE  FOR  COMPLEX  SOURCE  CALCULATIONS 

The  method  described  in  Appendix  A  provides  a  means  for 
computing  the  synthetic  seismograms  from  the  output  of  complex 

source  calculations.  Such  seismograms  include  both  source  and  path 
effects.  It  would  be  useful  to  isolate  the  source  effects  to  the 
extent  possible.  One  way  to  do  this  is  to  compare  two  seismograms 
computed  with  a  simple  source  and  the  same  path.  Differences 

between  the  two  are  then  due  to  differences  in  the  source. 

For  axisymmetric  explosion  calculations,  the  interesting 
comparison  is  to  a  spherical  point  source  representation  for  the 
explosion.  Let  the  source  function  for  such  an  explosion  be  the 
reduced  velocity  potential  (<u).  Then  the  Rayleigh  waves  for  this 
explosion  may  be  written  (e.g.,  Bache,  Rodi  and  Harkrider,  1978), 


WEX  * 


(B.l) 


where  P(r,u)  represents  the  path  effects.  Similarly,  the  Rayleigh 
waves  computed  with  the  wave  field  continuation  method  may  be  written 


WAC  "  'i'e  ’ 


(B.2) 


where  f  represents  some  "equivalent  RDP"  representation  for  the 
complex  source  calculation.  We  have 


As  an  example,  consider  the  test  problem  of  Figure  A. 3.  Them 
for  this  calculation  was  plotted  in  Figure  A. 2  and  ratio  w^^EX 
appears  in  Figure  A. 5.  In  Figure  B.l  we  show  the  Inland  |wgj  for 
this  comparison.  If  there  were  no  numerical  errors,  the  two  would 
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be  identical.  On  the  other  hand,  if  wp,(  had  been  computed  with  r 
at  a  different  depth  than  the  0.15  kilometer  depth  used  for 

A 

computing  the  monitoring  solutions,  we  would  expect  the  7  to  be 

A  ® 

much  less  close  to  ¥  . 

In  analyzing  nonlinear  axisymmetric  finite  difference 
calculations  of  explosions,  we  would  like  to  compare  the  wAC  to 
w^x  computed  with  a  ¥  from  a  one-dimensional  calculation  in  the 
same  material.  The  differences  between  the  two  are  then  precisely 
those  differences  caused  by  two-dimensional  effects. 
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